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I. INTRODUCTION 

The theoretical study of planet formation has a long 
history. Many of the fundamental ideas in the theory 
of terrestrial planet formation were laid out by Safronov 
(1969) in his classic monograph "Evolution of the Pro- 
toplanetary Cloud and Formation of the Earth and the 
Planets". The core accretion theory for gas giant for- 
mation, which was discussed by Cameron in the early 
1970's (Perri & Cameron, 1973), had been quantitatively 
developed in recognizable detail by 1980 (Mizuno, 1980). 
A wealth of new data over the last fifteen years — in- 
cluding direct observations of protoplanetary disks, the 
discovery of the Solar System's Kuiper Belt, and the de- 
tection of numerous extrasolar planetary systems — has 
led to renewed interest in the problem. Although these 
observations have confirmed some existing predictions, 
they have also emphasized the need to explore new the- 
oretical avenues. The major questions that work in this 
field seeks to answer include: 
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TABLE I Basic properties of planets in the Solar System, the 
semi- major axis a, eccentricity e, and mass M p . 





a/AU 


e 


M P /g 


Mercury 


0.387 


0.206 


3.3 x 10 2b 


Venus 


0.723 


0.007 


4.9 x 10 27 


Earth 


1.000 


0.017 


6.0 x 10 27 


Mars 


1.524 


0.093 


6.4 x 10 26 


Jupiter 


5.203 


0.048 


1.9 x 10 30 


Saturn 


9.537 


0.054 


5.7 x 10 29 


Uranus 


19.189 


0.047 


8.7 x 10 28 


Neptune 


30.070 


0.009 


1.0 x 10 29 



• How do terrestrial and giant planets form? 

• How much evolution in the orbits of planets takes 
place at early times? 

• Is the architecture of the Solar System typical? 

• How common are habitable planets, and are they 
inhabited? 

The goal of these notes is to provide a succinct introduc- 
tion to the concepts necessary to understand the astro- 
physics of planet formation. Before delving into theory, 
however, we first briefly review the basic observational 
properties of the Solar System and of extrasolar plan- 
etary systems that a theory of planet formation might 
aspire to explain. 

A. Critical Solar System Observations 

1. Architecture 

The orbital properties and masses of the planets in 
the Solar System are listed in Table I (the values here 
are taken from the JPL web site). The basic architec- 
ture of our Solar System comprises 2 gas giants. These 
are Jupiter and Saturn, which are composed primarily 
of hydrogen and helium - like the Sun - though they 
have a higher abundance of heavier elements as com- 
pared to Solar composition. Saturn is known to have 
a substantial core. Descending in mass there are then 
2 ice giants (Uranus and Neptune) composed of water, 
ammonia, methane, silicates and metals, plus low mass 
hydrogen / helium atmospheres; 2 large terrestrial plan- 
ets (Earth and Venus) plus two smaller terrestrial plan- 
ets (Mercury and Mars). Apart from Mercury, all of the 
planets have low eccentricities and orbital inclinations. 
They orbit in a plane that is approximately perpendic- 
ular to the Solar rotation axis (though there is a small 
misalignment of 7° ) . 

In the Solar System the giant and terrestrial planets 
are clearly segregated in orbital radius, with the inner 
zone occupied by the terrestrial planets being separated 
from the outer giant planet region by the main asteroid 
belt. The orbital radii of the giant planets coincide with 



where we expect the protoplanetary disk to have been 
cool enough for ices to have been present. This is a sig- 
nificant observation in the classical theory of giant planet 
formation, since in that theory the time scale for giant 
planet formation depends upon the mass of condensible 
materials. One would therefore expect faster growth to 
occur in the outer ice-rich part of the protoplanetary disk. 

2. Mass and angular momentum 

The mass of the Sun is M Q = 1.989 x 10 33 g, made 
up of hydrogen (fraction by mass X = 0.73), helium 
(Y = 0.25) and "metals" (Z = 0.02). One observes im- 
mediately that, 

M »^M p , (1) 

i.e. most of the heavy elements in the Solar System are 
found in the Sun rather than in the planets. The import 
of this trivial observation is that if most of the mass in 
the Sun passed through a disk during star formation, the 
planet formation process need not be very efficient. 

The angular momentum budget for the Solar System 
is dominated by the orbital angular momentum of the 
planets. The angular momentum in the Solar rotation is, 

L Q ~ k 2 M Q R%Q,, (2) 

assuming for simplicity solid body rotation. Taking 
O = 2.9 x 10~ 6 s- 1 and adopting k 2 = 0.1 (roughly 
appropriate for a star with a radiative core), L Q ~ 
3 x 10 48 gcm 2 s _1 . By comparison, the orbital angular 
momentum of Jupiter is, 

Lj = M Jy /GM Q a = 2 x lO^gcmV 1 . (3) 

The significance of this result is that it implies that sub- 
stantial segregation of mass and angular momentum must 
have taken place during (and subsequent to) the star for- 
mation process. We will look into how such segregation 
arises during disk accretion later. The broader question 
of exactly how the angular velocity of low mass stars 
evolves at early times remains a subject of active research 
(Herbst et al., 2007). 

3. Minimum mass Solar Nebula 

We can use the observed masses and compositions of 
the planets to derive a lower limit to the amount of gas 
that must have been present when the planets formed. 
This is called the Minimum Mass Solar Nebula (Weiden- 
schilling, 1977). The procedure is: 

1. Start from the known mass of heavy elements (say 
iron) in each planet, and augment this mass with 
enough hydrogen and helium to bring the mixture 
to Solar composition. This is a mild augmentation 
for Jupiter, but a lot more for the Earth. 
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2. Then divide the Solar System into annuli, with one 
planet per annulus. Distribute the augmented mass 
for each planet uniformly across the annuli, to yield 
a characteristic gas surface density £ (units gem -2 ) 
at the location of each planet. 

The result is that between Venus and Neptune (and 
ignoring the asteroid belt) S oc r~ 3 / 2 . To derive a pre- 
cise normalization from such a hand-waving procedure 
is somewhat pointless, but if one needs a specific num- 
ber the most common value used is that due to Hayashi 
(1981), 

s = L7x lo3 (xu)~ 3/2 gcm ~ 2 - W 

Integrating this expression out to 30 AU the enclosed 
mass works out to be around 0.01 M Q , comparable to the 
"typical" mass estimated for protoplanetary disks around 
other stars as inferred from mm observations of the dust. 

As the name should remind you, this is a minimum 
mass. It is not an estimate of the disk mass at the time 
the Solar Nebula formed, nor is there any reason to be- 
lieve that the S oc r~ 3 / 2 scaling represents the actual 
surface density profile for a protoplanetary disk. Most 
theoretical models of disks predict a significantly shal- 
lower slope more akin to E oc r _1 (Bell et al., 1997). 

4. Resonances 

A resonance occurs when there is a near-exact rela- 
tion between characteristic frequencies of two bodies. For 
example, a mean-motion resonance occurs between two 
planets with orbital periods Pi and Pi when, 




with i, j integers (the resonance is typically important if i 
and j are small integers). The "approximately equal to" 
sign in this expression reflects the fact that resonances 
have a finite width, which varies with the particular res- 
onance and with the eccentricities of the bodies involved. 
Resonant widths can be calculated accurately, but not 
easily (Murray & Dcrmott, 1999). In the Solar System 
Neptune and Pluto (along with many other Kuiper Belt 
objects) are in a 3:2 resonance, while Jupiter and Sat- 
urn are close to a 5:2 mean-motion resonance (known 
as "the great inequality") which influences their motion 
(Lovett, 1895). There are no simple resonances among 
the major planets. There are, however, many resonant 
pairs among planetary moons. Jupiter's satellites Io, Eu- 
ropa and Ganymede, for example, form a resonant chain 
in which Io is in 2:1 resonance with Europa, which itself 
is in a 2:1 resonance with Ganymede. The existence of 
these non-trivial configurations is generally regarded as 
strong circumstantial evidence that dissipative processes 
(such as tides) resulted in orbital evolution and trapping 
into resonance at some point in the past history of the 
systems (Goldreich, 1965). 



5. Minor bodies 

As a rough generalization the Solar System is dynam- 
ically full, in the sense that most locations where test 
particle orbits would be stable for 5 Gyr are, in fact, 
occupied by minor bodies. In the inner and middle So- 
lar System the main asteroid belt is the largest reservoir 
of minor bodies. The asteroid belt displays considerable 
structure, most notably in the form of sharp decreases in 
the number of asteroids in zones known as the Kirkwood 
gaps. The existence of these gaps provides a striking illus- 
tration of the importance of resonances (in this case with 
Jupiter) in influencing dynamics. Their detailed shapes 
provide constraints on the extent to which the orbits of 
the giant planets could have evolved during the lifetime 
of the Solar System (Minton & Malhotra, 2009). 

The properties of objects beyond Neptune (Chiang et 
al., 2007; Jewitt & Luu, 1993) provide further constraints 
on both the early evolution of the outer Solar System 
(Malhotra, 1993), and on collisional models for planet 
formation (Kenyon, 2002). Kuiper Belt properties in- 
clude: 

1 . A large population of objects in Pluto- like orbits in 
3:2 resonance with Neptune ( "plutinos" ) . 

2. A dearth of KBOs in orbits with 36 AU < a < 
39 AU. 

3. An apparent edge to the distribution of Classical 
KBOs at about 50 AU (Trujillo, Jewitt & Luu, 
2001). 

4. A differential size distribution (deduced indirectly 
from the measured luminosity function) that is 
roughly a power-law for large bodies with diame- 
ters D > 100 km (Trujillo, Jewitt & Luu, 2001). A 
recent determination by Fraser & Kavelaars (2009) 
infers a power-law slope q ~ 4.8 for large bodies 
together with a break to a much shallower slope at 
small sizes. 

Kuiper Belt Objects are commonly classified into several 
dynamically distinct families. Resonant KBOs are those 
- like Pluto — that exhibit mean-motion resonances 
with Neptune. Centaurs are non-resonant KBOs which 
have perihelion distances interior to the orbit of Neptune. 
Classical KBOs are objects further out whose orbits have 
been little influenced by Neptune. Finally, scattered disk 
KBOs are bodies with perihelia beyond the orbit of Nep- 
tune that do not fall into the other classes. 

Beyond the Kuiper belt the most intriguing object 
known is Sedna, a large object with semi-major axis 
a = 480 ± 40 AU, eccentricity e = 0.84 ± 0.01, and in- 
clination i = 12° (Brown, Trujillo & Rabinowitz, 2004). 
Since Sedna was discovered close to perihelion, it is highly 
likely to represent the first of a substantial new class of 
objects with perihelion distances substantially exterior to 
the orbit of Neptune. It may represent an object in an 
inner extension of the Oort cloud. 
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6. Ages 

Radioactive dating of meteorites provides an absolute 
measure of the age of the Solar System, together with 
constraints on the time scales of some phases of planet 
formation. The details are beyond the scope of these lec- 
tures; typical numbers quoted are a Solar System age of 
4.57 Gyr, a time scale for the formation of large bod- 
ies within the asteroid belt of < 5 Myr (Wadhwa et al., 
2007), and a time scale for final assembly of the Earth of 
~ 100 Myr. 



7. Satellites 

Most of the planets possess satellite systems, some of 
which are very extensive. Their origins appear to be 
heterogeneous. The regular satellites of the giant plan- 
ets have prograde orbits approximately coincident with 
the equatorial plane of their host planet. They appear 
to have formed in situ (though there may have been 
orbital evolution) from sub-disks of gas and dust that 
surrounded the planets at early times (Canup & Ward, 
2002, 2008). The irregular satellites 1 , on the other hand, 
appear to be captured objects (Jewitt & Haghighipour, 
2007). Finally, yet other satellites, most notably the 
Moon, accreted following giant impacts during planet for- 
mation (Benz, Slattery & Cameron, 1986; Canup, 2004). 



B. Extrasolar Planets 

1. Detection methods and biases 

The most important current methods for detecting and 
characterizing extrasolar planets are: 

1. Radial velocity surveys of nearby, typically Solar- 
type stars (Butler et al., 1996). More than 400 
planets have been found with this technique. 

2. Blind transit searches, and follow-up of radial ve- 
locity discovered planets that happen to show tran- 
sits (Charbonneau et al., 2007). There has been 
explosive growth in the success of this technique 
over the last couple of years, and more than 80 
such planets are now known. A number of ground- 
based surveys have now been joined by two space 
missions, COROT (Auvergne et al., 2009) and the 
larger and more sensitive Kepler (Borucki et al., 
2010). The realized photometric performance of 
Kepler should be sufficient to measure the stellar 



One conventional definition is that irregular satellites are those 
that orbit more than 0.05 Hill radii away from their planet, where 
the Hill radius is defined as r H = (M p /3M ) 1 / 3 a. 




FIG. 1 A planet of mass M p orbits the common center of mass 
at distance a\ , while the star of mass M* orbits at distance 
ct2. The system is observed at inclination angle i. 

flux decrement, 

caused by an Earth-radius planet orbiting a Solar- 
type star, though independent confirmation of the 
planetary origin of such a signal will be challeng- 
ing. Detection of small planets via transits from 
the ground is also possible, by focusing on low mass 
stars which have smaller physical radii (Charbon- 
neau et al., 2009). 

3. Gravitational lensing (Beaulieu et al., 2006), which 
is the ground-based technique with the greatest 
sensitivity for low mass planets. 

4. Direct imaging, which although presently limited to 
searches for massive planets at large orbital radii, 
has already uncovered the exceptionally interest- 
ing multiple system surrounding the star HR 8799 
(Marois et al., 2008). This system consists of three 
very massive planets orbiting at projected separa- 
tions between 20 and 70 AU from the star. 

5. Pulsar timing (Wolszczan & Frail, 1992). 

In addition to these methods astrometry (McArthur et 
al., 2010), and transit timing (Agol et al., 2005; Hol- 
man & Murray, 2005) have significant potential for both 
the discovery and characterization of extrasolar planetary 
systems. 

For statistical studies, the most important of the cur- 
rent methods is radial velocity surveys. 51 Peg, the first 
known extrasolar planet orbiting a normal star, was dis- 
covered this way (Mayor & Queloz, 1995), and most of 
our current knowledge of the extrasolar planet popula- 
tion derives from radial velocity surveys (Marcy et al., 
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FIG. 2 Highly schematic illustration of the selection function 
of an idealized radial velocity survey. The minimum mass 
planet that can be detected scales with semi-major axis as 
a 1 / 2 until the orbital period of the planet exceeds the duration 
of the survey. 



2005). The observable is the time dependence of the ra- 
dial velocity of a star due to the presence of an orbiting 
planet. For a planet on a circular orbit the geometry is 
shown in Figure 1. The star orbits the center of mass 
with a velocity, 



(7) 



Observing the system at an inclination angle i, we see the 
radial velocity vary with a semi-amplitude K = v* sini, 



K oc M p sin ia 



-1/2 



(8) 



If the inclination is unknown, what we measure (K) de- 
termines a lower limit to the planet mass M p . Note 
that is not determined from the radial velocity curve, 
but must instead be determined from the stellar spectral 
properties. If the planet has an eccentric orbit, e can be 
determined by fitting the non-sinusoidal radial velocity 
curve. 

The noise sources for radial velocity surveys comprise 
photon noise, intrinsic jitter in the star (e.g. from con- 
vection or stellar oscillations), and instrumental effects. 
The magnitude of these effects vary (sometimes dramat- 
ically) from star to star. However, if we imagine an ide- 
alized survey for which the noise per observation was a 
constant, then the selection limit would be defined by, 



Alp sin i | minimum 



Ca 1 / 2 , 



(9) 



with C a constant. Planets with masses below this 
threshold would be undetectable, as would planets with 
orbital periods exceeding the duration of the survey 
(since orbital solutions are poorly constrained when only 
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wavelength A 

FIG. 3 Schematic spectrum in the vicinity of a single spectral 
line of the host star. The wavelength range that corresponds 
to a single pixel in the observed spectrum is shown as the 
vertical shaded band. If the spectrum shifts by a velocity 5v 
the number of photons detected at that pixel will vary by an 
amount that depends upon the local slope of the spectrum. 



part of an orbit is observed unless the signal to noise of 
the observations is very high). The selection boundary 
defined by these limits is shown schematically in Figure 2. 

Extremely accurate radial velocity measurements are 
a prerequisite for discovering planets via this technique. 
For the Solar System, 



w 12 ms 1 (Jupiter) 
w 0.1 mr 1 (Earth). 



(10) 



Given that astronomical spectrographs have a resolving 
power of the order of 10 5 (which corresponds, in velocity 
units, to a precision of the order of kilometers per second) 
it might seem impossible to find planets with such small 
radial velocity signatures. To appreciate how detection of 
small (sub-pixel) shifts is possible, it is useful to consider 
the precision that is possible against the background of 
shot noise (i.e. uncertainty in the number of photons due 
purely to counting statistics). An estimate of the photon 
noise limit can be derived by considering a very simple 
problem: how accurately can velocity shifts be estimated 
given measurement of the flux in a single pixel on the de- 
tector? To do this, we follow the basic approach of Butler 
et al. (1996) and consider the spectrum in the vicinity of 
a spectral line, as shown in Figure 3. Assume that, in 
an observation of some given duration, -/V pn photons are 
detected in the wavelength interval corresponding to the 
shaded vertical band. If we now imagine displacing the 
spectrum by an amount (in velocity units) 5v the change 
in the mean number of photons is, 



SN ph 



dN ph 
dv 



Sv. 



(11) 



Since a la detection of the shift requires that 6N p h w 

1 /2 

iV pn , the minimum velocity displacement that is de- 
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tectable is, 



Sv n 



N. 



1/2 



ph 



dN ph /dv ■ 



(12) 



This formula makes intuitive sense - regions of the spec- 
trum that are fiat are useless for measuring Sv while sharp 
spectral features are good. For Solar-type stars with pho- 
tospheric temperatures T e g sa 6000 K the sound speed at 
the photosphere is around 10 kms -1 . Taking this as an 
estimate of the thermal broadening of spectral lines, the 
slope of the spectrum is at most, 



1 dK 



P h 



1 



N ph dv 



10 kms 



10~ 4 m _1 s. 



(13) 



Combining Equations (12) and (13) allows us to estimate 
the photon-limited radial velocity precision. For exam- 
ple, if the spectrum has a signal to noise ratio of 100 (and 
there are no other noise sources) then each pixel receives 
7V p h ~ 10 4 photons and Sv m i n ~ 100 ms _1 . If the spec- 
trum contains N p - lx such pixels the combined limit to the 
radial velocity precision is, 



5v< 



8v n 



100 ms" 



hot 



N 



1/2 



N. 



1/2 



(14) 



Obviously this discussion ignores many aspects that are 
practically important in searching for planets from ra- 
dial velocity data. However, it suffices to reveal the key 
feature: given a high signal to noise spectrum and sta- 
ble wavelength calibration, photon noise is small enough 
that a radial velocity measurement with the ms _1 preci- 
sion needed to detect extrasolar planets is feasible. 

Records for the smallest amplitude radial velocity sig- 
nal that can be extracted from the noise are regularly 
bested. Currently some of the highest precision radial 
velocity measurements have an RMS scatter of around 
0.5 ms _1 , and it has been suggested that 0.2 ms _1 might 
be attainable for some quiescent stars (Mayor & Udry, 
2008). The lowest stellar velocity semi- amplitude is 
slightly below 2 ms" 1 . It is important to remember that 
these are best-case values - complete samples of extra- 
solar planets that are suitable for statistical studies only 
exist for much larger K 30 ms _1 (Fischer & Valenti, 
2005). 

Detailed modeling is necessary in order to assess 
whether a particular survey has a selection bias in eccen- 
tricity. Naively you can argue it either way - an eccentric 
planet produces a larger perturbation at closest stellar 
approach, but most of the time the planet is further out 
and the radial velocity is smaller. A good starting point 
for studying these issues is the explicit calculation for the 
Keck Planet Search reported by Cumming et al. (2008). 
These authors find that the Keck search is complete for 
sufficiently massive planets (and thus trivially unbiased) 
for e < 0.6. 
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FIG. 4 The distribution of known extrasolar planets in semi- 
major axis and eccentricity (red triangles). Solar System 
planets are shown for comparison as the blue squares. The 
dashed curve denotes a line of constant periastron distance. 
The figure uses data from an updated version of the But- 
ler et al. (2006) catalog, and includes planets that have 
Mpsini < lOMj. 



2. Observed properties 

For most known extrasolar planets, our information is 
limited to those quantities derived from the radial veloc- 
ity observables: a lower limit on the mass M p sini, the 
semi-major axis a, the eccentricity e, and the longitude 
of pericenter w. In addition, estimates of the host star's 
mass and metallicity are available. The distribution of 
planets in M p sin i, a and e is depicted in Figures 4, 5 and 
6, using data for radial velocity detected planets from an 
updated version of the Butler et al. (2006) catalog. 

Marcy et al. (2005) quote the following results from 
the Lick / Keck / AAT survey, which has monitored 1330 
FGKM stars for the better part of a decade: 

1. The detected giant planet frequency within a ~ 
5 AU is ~ 7%. This is certainly a lower limit as 
many giant planets would fall below the selection 
threshold at larger orbital radii (c.f. Figure 5). 

2. The frequency of hot Jupiters with a < 0.1 AU is 
approximately 1%. The abundance of planets with 



orbital radius 
to large a. 



measured as dN p /d log a - increases 



3. Eccentric orbits are common beyond the radius 
where tidal circularization is significant (Figure 4) . 
The median eccentricity of planets orbiting between 
1 and 3 AU is (e) ~ 0.28. Some extremely eccentric 
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FIG. 5 The distribution of known extrasolar planets, discov- 
ered via various techniques, is shown in the space of semi- 
major axis and minimum mass. A line of constant semi- 
amplitude radial velocity perturbation (K = 5 m s" 1 ) is plot- 
ted assuming a Solar mass host. It is clear by eye that the 
typical extrasolar planet detected so far is not a hot Jupiter 
but rather orbits at a > 1 AU. 



1 - 



0.8 - 



0.6 



o 0.4 

CD 



0.2 - 



Short period planets (a < 0.1 AU) 
All other systems ■ 



■ \ ■ 



•- 



. . ■ ■ _ ■ ■ 

" * . ■-■ ■ . " 



0.1 1 10 

M sin (i) / Jupiter masses 



FIG. 6 Eccentricity vs mass for known extrasolar planets, di- 
vided into short period planets (a < 0.1 AU, shown as blue 
triangles) and all other systems (shown as red squares). The 
short period planets have low eccentricity orbits, presumably 
as a consequence of tidal interactions with their host stars. 
Although the eccentricity distribution is not entirely inde- 
pendent of mass the correlation is weak. 
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FIG. 7 The fraction of stars that host currently known extra- 
solar planets is plotted as a function of the stellar metallic- 
ity, from data (their Figure 4) reported by Fischer & Valenti 
(2005). 



planets exist. There is no strong trend of eccentric- 
ity with planet mass (Figure 6). 

4. The planet mass function declines toward large 
masses (Butler et al., 2006; Tabachnik & Tremaine, 
2002). 

5. Planet frequency rises rapidly with host metallicity. 
This trend, shown in Figure 7 using data from Fis- 
cher & Valenti (2005) is dramatic — relatively mod- 
est increases in metallicity substantially enhance 
the probability that currently detectable planets 
will be found around a star. 

6. Multiple planet systems are common, of which a 
significant fraction exhibit prominent mean-motion 
resonances between planets. 

Transit surveys contribute complementary information 
about the extrasolar planet population. Ground based 
surveys - which due to atmospheric limitations on pho- 
tometric precision are only sensitive to massive plan- 
ets - have discovered numerous short and ultrashort 
period planets (e.g. WASP-12b with a period of only 
1.09 days) and measured their physical radii. The plan- 
etary radii confirm their gas giant nature. Space based 
transit surveys, together with the highest precision radial 
velocity measurements, have extended this sample into 
the so-called super-Earth regime of planets with masses 
M p < 10 M®. The properties of the super-Earth popu- 
lation remain somewhat murky, but these may plausibly 
be rocky (or in one case water-rich) bodies that bridge 
the gap between the Solar System's terrestrial and giant 
planets. 

Transit observations have also provided two surprises. 
First, although the basic confirmation that massive plan- 
ets have radii similar to that of Jupiter is reassuring, 
the scatter in the observed radii is inconsistent with the 
simplest theoretical predictions. These radius anoma- 
lies (Torres, Winn & Holman, 2008) comes in two flavors 
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(Burrows et al., 2007): some planets are "too small" (as 
compared to theoretical models), while others are "too 
large". The undersized gas giants are interesting, but 
pose no special theoretical conundrum. To first order, 
the radius of a gas giant of a given mass varies with the 
total mass of heavy elements it contains 2 ; hence a plau- 
sible explanation for any small planet is that it has an 
above-average heavy element content. The measured ra- 
dius of the Saturn mass planet orbiting HD 149026, for 
example, is generally interpreted as providing evidence 
for approximately 70M e of heavy elements in the inte- 
rior (Sato et al., 2005). The inflated planets, on the other 
hand, are more mysterious, since some (most particularly 
TrES-4 and WASP- 12b) are too large even when com- 
pared to pure hydrogen / helium models. It appears as 
if explaining their radii requires an additional source of 
heat, which could in principle come from several sources, 

• Stellar irradiation. This is the obvious answer, 
since the planets with inflated radii all have very 
short-period orbits, and only a small fraction of 
the stellar insolation would suffice to explain the 
measured radii. The difficulty is that the irradiat- 
ing flux only directly heats the surface, while the 
radius depends upon the entropy of the convective 
core that does not extend to the surface. A direct 
solution that appeals to irradiation is only viable if 
the opacity in the atmosphere is substantially dif- 
ferent from that predicted by current models (Bur- 
rows et al., 2007). 

• Atmosphere-interior coupling. Since the magnitude 
(if not the location) of stellar irradiation is more 
than adequate, we can tolerate some inefficiency 
and postulate a mechanism that mixes some frac- 
tion of the energy from the surface into the convec- 
tive zone. Guillot & Showman (2002), for example, 
suggested that some of the energy in surface winds 
(driven by irradiation) could mix into the interior 
via waves. More recently, Batygin & Stevenson 
(2010) have shown that plausible magnetic fields 
could furnish the required coupling. In this model, 
the interaction of surface winds with the magnetic 
field generates currents in the planetary interior, 
which are dissipated resistively. 

• Tides. A planet with a non-zero eccentricity will be 
heated by tides raised on the planet by the nearby 
star. This mechanism must be important for some 
or all planets - particularly at early times - but 
there are reasons to doubt that it is the sole expla- 
nation for inflated planetary radii. In particular, 
there are stringent limits in some systems to the 
presence of additional planets that could maintain 



2 Whether those heavy elements are distributed evenly within the 
planet, or concentrated at the center in a core, also affects the 
radius, but at a more subtle level. 



a finite eccentricity against the dissipativc effects 
of the tidal interactions themselves. 

Observationally, it may be possible to divine which of 
these mechanisms (or, some other possibility) is at work 
by searching for correlations between radius anomalies 
and other planetary properties - such as the existence of 
temperature inversions in the atmosphere or additional 
planets in the same system. 

In addition to the longstanding mystery of the radii, 
transits have recently uncovered a second puzzle. In a 
number of systems the combination of transit and ra- 
dial velocity data has permitted a measurement of the 
orientation of the orbital plane relative to the stellar 
spin axis via detection of the Rossiter-McLaughlin effect 
(McLaughlin, 1924; Rossiter, 1924) 3 . Both intuition and 
Solar System precedent would suggest that any misalign- 
ment between the stellar and orbital angular momentum 
vectors ought to be small, but while this is often the case 
at least some systems show quite large misalignments. 
The XO-3 system, for example, shows a well measured 
misalignment of 37±4 degrees (Winn et al., 2009). This 
is of interest for the possible constraints it places on the- 
oretical models of orbital migration of massive planets. 



II. PROTOPLANETARY DISKS 

A. The star formation context 

Stars form in the Galaxy today from the small fraction 
of gas that exists in dense, molecular clouds. Molecular 
clouds are observed in one or more molecular tracers - 
examples include CO, 13 CO and NH3 - which can be 
used both to probe different regimes of column density 
and to furnish kinematic information that can give clues 
as to the presence of rotation, infall and outflows. Obser- 
vations of the dense, small scale cores within molecular 
clouds (with scales of the order of 0.1 pc) that are the 
immediate precursors of star formation show velocity gra- 
dients that are of the order of 1 kms _1 pc _1 . Even if all of 
such a gradient is attributed to rotation, the parameter, 



|-^grav| 

is small - often of the order of 0.01. Hence rotation is 
dynamically unimportant during the early stages of col- 
lapse. The angular momentum, on the other hand, is 



The effect consists of a shift in the apparent stellar radial ve- 
locity during transit as the planet obscures parts of the pho- 
tosphere that arc rotating cither toward or away from the ob- 
server. Clear detections of the "rotational effect" (as it was then 
called) in eclipsing binaries were published by Richard Rossiter 
(as part of his Ph.D. studying the beta Lyrae system), and by 
Dean McLaughlin (who studied Algol). Frank Schlesinger, and 
possibly others, may have seen similar effects in binaries. 
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,1 micron 




to define a measure of the slope of the IR SED, 



log X 

FIG. 8 Schematic depiction of the Spectral Energy Distribu- 
tion of a young star surrounded by a disk. The presence of a 
disk is inferred from an infra-red excess (above the expected 
photospheric value) at wavelengths longward of around 1 /im. 
An ultra-violet excess is also commonly detected, and this is 
attributed to gas accretion on to the stellar surface producing 
hot spots. 



large, with a ballpark figure being J core ~ 10 54 gcm 2 s _1 . 
This is much larger than the angular momentum in the 
Solar System, never mind that of the Sun, a discrepancy 
that is described as the angular momentum problem of 
star formation. The overall solution to this problem is 
thought to be an undetermined admixture of binary for- 
mation, angular momentum loss in outflows, and disk for- 
mation. For our purposes, it suffices to note that the spe- 
cific angular momentum of gas in molecular cloud cores 
would typically match the specific angular momentum 
of gas in Keplerian orbit around a Solar mass star at a 
radius of - 10 - 10 2 AU. 

The bottom line is thus simply that the observed prop- 
erties of molecular cloud cores are consistent with the 
formation of large disks - of the size of the Solar System 
and above - around newly formed stars. At least initially, 
those disks could be quite massive. 

Young Stellar Objects (YSOs) are classified observa- 
tionally according to the shape of their Spectral Energy 
Distribution XF\(X) in the infra-red. As shown schemat- 
ically in Figure 8, YSOs often display, 

1. An infra-red excess (over the stellar photospheric 
contribution) - this is attributed to hot dust in the 
disk near the star. 

2. An ultra-violet excess, which is ascribed to high 
temperature regions (probably hot spots) on the 
stellar surface where gas from the disk is being ac- 
creted. 

To quantify the magnitude of the IR excess, it is useful 



am 



Alog(Af\) 
A log A 



(16) 



between the near-IR and the mid-IR. Conventions vary, 
but for illustration we can assume that the slope is mea- 
sured between the K band (at 2.2/im) and the N band 
(at 10/zm). We can then classify YSOs as, 

• Class 0: SED peaks in the far-IR or mm part of 
the spectrum (~ 100 /Ma), with no flux being de- 
tectable in the near-IR. 

• Class I: approximately flat or rising SED into mid- 
IR (am > 0). 

• Class II: falling SED into mid-IR (—1.5 < axR < 
0). These objects are also called "Classical T Tauri 
stars" . 

• Class III: pre-main-sequence stars with little or 
no excess in the IR. These are the "Weak lined 
T Tauri stars" (note that although WTTs are de- 
fined via the equivalent width of the Ha line, this 
is an accretion signature that correlates well with 
the presence of an IR excess). 

This observational classification scheme is theoretically 
interpreted, in part, as an evolutionary sequence (Adams, 
Lada & Shu, 1987). In particular, clearly objects in 
Classes through II eventually lose their disks and be- 
come Class III sources. Viewing angle may well, however, 
play a role in determining whether a given source is ob- 
served as a Class I or Class II object. 



B. Passive circumstellar disks 

An important physical distinction needs to be drawn 
between passive circumstellar disks, which derive most 
of their luminosity from reprocessed starlight, and active 
disks, which are instead powered by the release of gravi- 
tational potential energy as gas flows inward. For a disk 
with an accretion rate M, surrounding a star with lumi- 
nosity Lq and radius i?* = 2i? Q , the critical accretion 
rate below which the accretion energy can be neglected 
may be estimated as, 



1 

4 L ° 



GNhM 
2i?„ 



(17) 



where we have anticipated the result, derived below, that 
a flat disk intercepts one quarter of the stellar flux. Nu- 
merically, 



M 



3 x 10" 8 M yr _x . 



(18) 



Measured accretion rates of Classical T Tauri stars (Gull- 
bring et al., 1998) range from an order of magnitude 
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FIG. 9 Geometry for calculation of the vertical hydrostatic 
equilibrium of a circumstellar disk. 



above this critical rate to two orders of magnitude be- 
low, so it is oversimplifying to assume that protoplan- 
etary disks are either always passive or always active. 
Rather, the dominant source of energy for a disk is likely 
to be a function of both time and radius. We expect 
internal heating to dominate at early epochs and / or 
small orbital radii, while at late times and at large radii 
reprocessing dominates. 



1. Vertical structure 

The vertical structure of a geometrically thin disk (ei- 
ther passive or active) is derived by considering vertical 
hydrostatic equilibrium (Figure 9). The pressure gradi- 
ent, 



dP 

"d7 



= ~P9z 



(19) 



where p is the gas density. Ignoring any contribution 
to the gravitational force from the disk (this is justified 
provided that the disk is not too massive), the vertical 
component of gravity seen by a parcel of gas at cylindrical 
radius r and height above the midplane z is, 



GM , 
d 2 



■ sin 6 



GM, 
d 3 



For a thin disk z <C r, so 



n 2 



(20) 



(21) 



where ft = y/ GM, / r 3 is the Keplerian angular veloc- 
ity. If we assume for simplicity that the disk is vertically 
isothermal (this will be a decent approximation for a pas- 
sive disk, less so for an active disk) then the equation of 
state is P = pc 2 , where c s is the (constant) sound speed. 
The equation of hydrostatic equilibrium (equation 19) 
then becomes, 



2 d P n 2 

c.^ = -(l pz. 



The solution is, 



2 /2h 2 



where h, the vertical scale height, is given by, 



(22) 



(23) 



(24) 



FIG. 10 Geometry for calculating the temperature profile of 
a flat, passive disk. We consider unit surface area in the disk 
plane at distance r from a star of radius f?». The axis of 
spherical polar co-ordinates is the line between the surface 
and the center of the star, with <j> — in the direction of the 
stellar pole. 



Comparing the thickness to the radius, 

- — — 

r W0 



(25) 



where v$ is the local orbital velocity. We see that the 
aspect ratio of the disk h/r is inversely proportional to 
the Mach number of the flow. 

The shape of the disk depends upon h(r)/r. If we pa- 
rameterize the radial variation of the sound speed via, 



c s oc r 

then the aspect ratio varies as, 
h 



-P 



/8+1/2 



(26) 



(27) 



The disk will flare - i.e. h/r will increase with radius 
giving the disk a bowl-like shape - if j3 < 1/2. This 
requires a temperature profile T(r) oc r _1 or shallower. 
As we will show shortly, flaring disks are expected to be 
the norm, at least relatively close to the star. 



2. Radial temperature profile 

The physics of the calculation of the radial tempera- 
ture profile of a passive disk is described in papers by 
Adams & Shu (1986), Kenyon & Hartmann (1987) and 
Chiang & Goldreich (1997). We begin by considering the 
absolute simplest model: a flat thin disk in the equato- 
rial plane that absorbs all incident stellar radiation and 
re-emits it as a single temperature blackbody. The back- 
warming of the star by the disk is neglected. 

We consider a surface in the plane of the disk at dis- 
tance r from a star of radius R, . The star is assumed to 
be a sphere of constant brightness /* . Setting up spheri- 
cal polar co-ordinates, as shown in Figure 10, the stellar 
flux passing through this surface is, 



F= I, sin 9 cos 4>dQ. 



(28) 



We count the flux coming from the top half of the star 
only (and to be consistent equate that to radiation from 
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only the top surface of the disk), so the limits on the 
integral are, 

-7I-/2 < 4> < tt/2 

< 6 < sin" 1 (^j . (29) 
Substituting dfl = sm8d9dcf>, the integral for the flux is, 

i/2 / .sin" 1 (i?.»/r) 



F = h , 

-7T/2 

which evaluates to, 



F = L 




(30) 



(31) 



For a star with effective temperature T„, the brightness 
/„ = (l/ir)aT 4 , with cr the Stefan-Boltzmann constant 
(Rybicki & Lightman, 1979). Equating F to the one- 
sided disk emission &T^ isk we obtain a radial temperature 
profile, 



?disk 



R* 
r 




Integrating over radii, we obtain the total disk flux, 



(32) 



F, 



disk 



2 x 



If, 

4 



2nraT disk dr 



(33) 



We conclude that a flat passive disk extending all the way 
to the stellar equator intercepts a quarter of the stellar 
flux. The ratio of the observed bolometric luminosity 
of such a disk to the stellar luminosity will vary with 
viewing angle, but clearly a flat passive disk is predicted 
to be less luminous than the star. 

The form of the temperature profile given by equation 
(32) is not very transparent. Expanding the right hand 
side in a Taylor series, assuming that (R*/r) <C 1 (i.e. 
far from the stellar surface), we obtain, 



Tdisk oc r 



-3/4 



(34) 



as the limiting temperature profile of a thin, flat, passive 
disk. For fixed molecular weight fi this in turn implies a 
sound speed profile, 

-3/8 



c s oc r 



(35) 

Assuming vertical isothermality, the aspect ratio given 
by equation (27) is, 



h 

— oc r 
r 



1/8 



(36) 



and we predict that the disk ought to flare modestly to 
larger radii. If the disk does flare, then the outer regions 
intercept a larger fraction of stellar photons, leading to 
a higher temperature. As a consequence, a temperature 
profile Idisk oc r~ 3 / 4 is probably the steepest profile we 
would expect to obtain for a passive disk. 



flat 'disk' 
part of SED 




log X 

FIG. 11 Schematic disk spectrum. At short wavelengths, we 
see an exponential cut-off corresponding to the highest tem- 
perature annulus in the disk (normally close to or at the inner 
edge). At long wavelengths, there is a Rayleigh- Jeans tail re- 
flecting the coldest material in the outer disk. At intermediate 
wavelengths, there is a flatter portion of the spectrum, so that 
the overall SED resembles a stretched blackbody. 



3. Spectral energy distribution (SED) 

Suppose that each annulus in the disk radiates as a 
blackbody at the local temperature IdiskM- If the disk 
extends from r m to r out , the disk spectrum is just the 
sum of these blackbodies weighted by the disk area, 



F x <x 



2irrB x [T(r)]dr 



where B\ is the Planck function, 



Bx(T) = 



2hc 2 



1 



\5 e hc/XkT _ I ' 



(37) 



(38) 



The behavior of the spectrum implied by equation (37) 
is easy to derive. At long wavelengths A ^> hc/kT(r out ) 
we recover the Rayleigh- Jeans form, 



XF X cx A~ 3 



(39) 



while at short wavelengths A <§; hc/kT(r- m ) there is an 
exponential cut-off that matches that of the hottest an- 
nulus in the disk, 



XFxot X - 4 e- hc / XkT(r -l 
For intermediate wavelengths, 



he 



«A« 



he 



(40) 



(41) 



kT(r in ) ' kT(r out ) 

the form of the spectrum can be found by substituting, 



he 



XkT(r b 



3/4 



(42) 
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into equation (37). We then have, approximately, 

x 5 / 3 d. 



and so 



Fx oc A- 7 / 3 ^ X -^- oc A- 7 /3 (43) 
AF A oc A" 4 / 3 . 



(44) 



The overall spectrum, shown schematically in Figure 11, 
is that of a "stretched" blackbody (Lynden-Bell, 1969). 

The SED predicted by this simple model generates an 
IR-excess, but with a declining SED in the mid-IR. This 
is too steep to match the observations of even most Class 
II sources. 



4. Sketch of more complete models 

Two additional pieces of physics need to be included 
when computing detailed models of the SEDs of passive 
disks. First, as already noted above, all reasonable disk 
models flare toward large r, and as a consequence inter- 
cept and reprocess a larger fraction of the stellar flux. At 
large radii, Kenyon & Hartmann (1987) find that consis- 
tent flared disk models approach a temperature profile, 



T disk cx r" 1 / 2 , 



(45) 



which is much flatter than the profile derived previously. 
Second, the assumption that the emission from the disk 
can be approximated as a single blackbody is too simple. 
In fact, dust in the surface layers of the disk radiates at a 
significantly higher temperature because the dust is more 
efficient at absorbing short-wavelength stellar radiation 
than it is at emitting in the IR (Shlosman & Begelman, 
1989). Dust particles of size a absorb radiation efficiently 
for A < 2ira, but are inefficient absorbers and emitters for 
A > 2na (i.e. the opacity is a declining function of wave- 
length). As a result, the disk absorbs stellar radiation 
close to the surface (where t ~ 1), where the optical 
depth to emission at longer IR wavelengths tir <C 1. The 
surface emission comes from low optical depth, and is not 
at the blackbody temperature previously derived. Chi- 
ang & Goldreich (1997) showed that a relatively simple 
disk model made up of, 

1. A hot surface dust layer that directly re-radiates 
half of the stellar flux 

2. A cooler disk interior that reprocesses the other 
half of the stellar flux and re-emits it as thermal 
radiation 

can, when combined with a flaring geometry, reproduce 
most SEDs quite well. A review of recent disk modeling 
work is given by Dullemond et al. (2007). 

The above considerations are largely sufficient to un- 
derstand the structure and SEDs of Class II sources. For 
Class I sources, however, the possible presence of an en- 
velope (usually envisaged to comprise dust and gas that 



is still infalling toward the star-disk system) also needs 
to be considered. The reader is directed to Eisner et al. 
(2005) for one example of how modeling of such systems 
can be used to try and constrain their physical properties 
and evolutionary state. 



C. Actively accreting disks 

The radial force balance in a passive disk includes con- 
tributions from gravity, centrifugal force, and radial pres- 
sure gradients. The equation reads, 



r 



GM, 



ldP 

p dr ' 



(46) 



where is the orbital velocity of the gas and P is the 
pressure. To estimate the magnitude of the pressure gra- 
dient term we note that, 



1 dP 

p dr 



IP 

p r 

p r 
GM* 



(47) 



where for the final step we have made use of the relation 
h = c 3 /ft. If vk is the Keplerian velocity at radius r, we 
then have that, 



2 2 

vs =v K 



I - 01* 

r 



(48) 



i.e pressure gradients make a negligible contribution to 
the rotation curve of gas in a geometrically thin (h/r <C 
1) disk 4 . To a good approximation, the specific angular 
momentum of the gas within the disk is just that of a 
Keplerian orbit, 



I = r 2 n = \J GM„r, 



(49) 



which is an increasing function of radius. To accrete 
on to the star, gas in a disk must lose angular momentum, 
either, 

1. Via redistribution of angular momentum within the 
disk (normally described as being due to "viscos- 
ity", though this is a loaded term, best avoided 
where possible). 

2. Via loss of angular momentum from the star-disk 
system, for example in a magnetically driven disk 
wind. 



4 This is not to say that pressure gradients arc unimportant - as 
we will see later the small difference between v$ and vk is of 
critical importance for the dynamics of small rocks within the 
disk. 
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Aspects of models in the second class have been stud- 
ied for a long time - the famous disk wind solution of 
Blandford & Payne (1982), for example, describes how a 
wind can carry away angular momentum from an under- 
lying disk. Observationally, it is not known whether mag- 
netic winds are launched from protoplanetary disks on 
1 — 100 AU scales (jets, of course, are observed, but these 
are probably launched closer to the star), and hence the 
question of whether winds are important for the large- 
scale evolution of disks remains open. An excellent re- 
view of the theory of disk winds, as applied to proto- 
stellar systems, is given by Konigl & Salmcron (2010). 
Here, we will assume that winds are not the dominant 
driver of evolution, and instead derive the equation for 
the time evolution of the surface density for a thin, vis- 
cous disk (Lynden-Bell & Pringle, 1974; Shakura & Sun- 
yaev, 1973). Clear reviews of the fundamentals of accre- 
tion disk theory can be found in Pringle (1981) and in 
Frank, King & Raine (2002). 



1. Diffusive evolution equation 

Let the disk have surface density S(r, t) and radial ve- 
locity v r (r,t) (defined such that v r < for inflow). The 
potential is assumed fixed so that the angular velocity 
ft = Q(r) only. In cylindrical co-ordinates, the conti- 
nuity equation for an axisymmetric flow gives (see e.g. 
Pringle (1981) for an elementary derivation), 

<9E d 

r ^ + d~r (r&V) = °- (50) 
Similarly, conservation of angular momentum yields, 



<9(Er 



7T ( rSlv 
or 



(51) 



dt dr 27r dr 

where the term on the right-hand side represents the net 
torque acting on the fluid due to viscous stresses. From 
fluid dynamics (Pringle, 1981), G is given in terms of the 
kinematic viscosity v by the expression, 



an 

G = 2rrr ■ z/Er— — ■ r 
dr 



(52) 



where the right-hand side is the product of the circum- 
ference, the viscous force per unit length, and the level 
arm r. If we substitute for G, eliminate v r between equa- 
tion (50) and equation (51), and specialize to a Keplerian 
potential with oc r~ 3 / 2 , we obtain the evolution equa- 
tion for the surface density of a thin accretion disk in its 
normal form, 



as 

~dt 



3d_ 

r dr 



■4 k /2 ) 



(53) 



This partial differential equation for the evolution of the 
surface density E has the form of a diffusion equation. 
To make that explicit, we change variables to, 

X = 2r^ 2 

f = ^X. (54) 



For a constant v, equation (53) then takes the prototyp- 
ical form for a diffusion equation, 



df _ a 2 f 
dt dx 2 



with a diffusion coefficient, 



D = 



\2v 

X2' 



(55) 



(56) 



The characteristic diffusion time scale implied by equa- 
tion (55) is X 2 1 D. Converting back to the physical vari- 
ables, we find that the evolution time scale for a disk of 
scale r with kinematic viscosity v is, 



(57) 



Observations of disk evolution (for example determina- 
tions of the time scale for the secular decline in the ac- 
cretion rate) can therefore be combined with estimates of 
the disk size to yield an estimate of the effective viscosity 
in the disk (Hartmann et al., 1998). 



2. Solutions 

In general, v is expected to be some function of the 
local conditions within the disk (surface density, radius, 
temperature, ionization fraction etc). If v depends on E, 
then equation (53) becomes a non-linear equation with 
no analytic solution (except in some special cases), while 
if there is more a complex dependence on the local con- 
ditions then the surface density evolution equation will 
often need to be solved simultaneously with an evolution 
equation for the central temperature (Pringle, Vcrbunt & 
Wade, 1986). Analytic solutions are possible, however, if 
v can be written as a power-law in radius (Lynden-Bell & 
Pringle, 1974), and these suffice to illustrate the essential 
behavior implied by equation (53). 

First, we describe a Green's function solution to equa- 
tion (53) for the case v = constant. Suppose that at 
t = 0, all of the gas lies in a thin ring of mass m at 
radius r , 



E(r,t = 0) 



27T7"0 



S(r - r ). 



One can show that the solution is then, 



X(x,r) = ~x- 1 / 4 e-( 1+ * 2 ^ Il 



/4 



(58) 



(59) 



12vr Q z t, and is a 



where we have written the solution in terms of dimen- 
sionless variables x = r/r a , r ~ 1 
modified Bessel function of the first kind 

Unless you have a special affinity for Bessel functions, 
this Green's function solution is not terribly transparent. 
The evolution it implies is shown in Figure 12. The most 
important features of the solution are that, as t — > oo, 
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FIG. 12 The Green's function solution to the disk evolution 
equation with v — constant, showing the spreading of a ring 
of mass initially orbiting at r = ro. From top down the curves 
show the behavior as a function of the scaled time variable 



t = 0.064, r 



for r = 0.004, r = 0.008, 
= 0.128, and r = 0.256. 



0.016, r = 0.032, 




100 



r / r, 



FIG. 13 The self-similar solution to the disk evolution equa- 
tion is plotted for a viscosity veer. The initial surface density 
tracks the profile for a steady-state disk (with E oc r~ L ) at 
small radius, before cutting off exponentially beyond r = ri. 
The curves show the surface density at the initial value of the 
scaled time T = 1, and at subsequent times T = 2, T = 4 and 
T = 8. 



dQ. /dr = 




stellar 
surface 

FIG. 14 Schematic depiction of the angular velocity Q(r) for 
a slowly rotating star surrounded by a thin accretion disk 
that extends to the stellar equator. At large radii in the disk, 
the angular velocity has the normal Keplerian form Q~ 3 ^ 2 , 
shown as the dashed green curve. To match smoothly on to 
the star, the angular velocity must turn over at smaller radii in 
a transition zone known as the boundary layer. The existence 
of a boundary layer implies that at some radius dSl/dr = 0, 
at which point the viscous stress vanishes. 



Suppose that the disk at time t = has the surface den- 
sity profile corresponding to a steady-state solution (with 
this viscosity law) out to r = r%, with an exponential cut- 
off at larger radii. As we will shortly show, the initial 
surface density then has the form, 



E(t = 0) 



C 



— exp 



_f(2-7) 



(61) 



where C is a normalization constant, f = r/r±, and i>\ 
v(ri). The self-similar solution is then, 



£(f,T) 



where, 



C 



37T^lf 7 



_ r -(5/2- 7 )/(2- 7 ) 



exp 



f(2-7) 



T 



(62) 



• The mass flows to r = 0. 

• The angular momentum, carried by a negligible 
fraction of the mass, flows toward r = oo. 

This segregation of mass and angular momentum is a 
generic feature of viscous disk evolution, and is obviously 
relevant to the angular momentum problem of star for- 
mation. 

Of greater practical utility is the self-similar solution 
also derived by Lynden-Bell & Pringle (1974). Consider 
a disk in which the viscosity can be approximated as a 
power-law in radius, 



T = 



v oc r ' 



(60) 



E + 1 



3(2 - 7 ) 2 v x 



(63) 



This solution is plotted in Figure 13. Over time, the disk 
mass decreases while the characteristic scale of the disk 
(initially r%) expands to conserve angular momentum. 
This solution is quite useful both for studying evolving 
disks analytically, and for comparing observations of disk 
masses, accretion rates or radii with theory (Hartmann 
et ah, 1998). 

A steady-state solution for the radial dependence of 
the surface density can be derived by setting d/dt = 
and integrating the angular momentum conservation 
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equation (51). This yields, 



,dfl 



'Er 3 Qv r = i^£r 3 — h constant. 

dr 



(64) 



Noting that the mass accretion rate M = — 2itrY 1 v r we 
have, 

-— r 2 n = ^Er 3 — + constant. (65) 
2tt dr v ' 

To determine the constant of integration, we note that 
the torque within the disk vanishes if dO/dr = 0. At 
such a location, the constant can be evaluated and is 
just proportional to the local flux of angular momentum 



constant oc Mr 2 Ct. 



(66) 



Usually this is determined at the inner boundary. A par- 
ticularly simple example is the case of a disk that extends 
to the equator of a slowly rotating star. This case is il- 
lustrated in Figure 14. In order for there to be a transi- 
tion between the Keplerian angular velocity profile in the 
disk and the much smaller angular velocity at the stellar 
surface there must be a maximum in it at some radius 
r* + Ar. Elementary arguments (Pringle, 1977) - which 
may fail at the very high accretion rates of FU Orionis 
objects (Popham ct al., 1993) but which are probably 
reliable otherwise - suggest that Ar <C r*, so that the 
transition occurs in a narrow boundary layer close to the 
stellar surface 5 . The constant can then be evaluated as, 



M 2 GM» 
constant*-— r, J —g- 



and equation (65) becomes, 



37T V V T 



(67) 



(68) 



Given a viscosity, this equation defines the steady-state 
surface density profile for a disk with an accretion rate 
M. Away from the boundaries, E(r) oc v^ 1 . 

The inner boundary condition which leads to equation 
(68) is described as a zero-torque boundary condition. 
As noted, zero-torque conditions are physically realized 
in the case where there is a boundary layer between the 
star and its disk. This is not, however, the case in most 
Classical T Tauri stars. Observational evidence suggests 
(Bouvier et al., 2007) that in accreting T Tauri stars 



5 The physics of the boundary layer itself presents interesting com- 
plications, since the boundary layer is a region of strong shear 
that is stable against the magnetorotational instabilities that we 
will argue later are critical for transporting angular momentum 
within disks. Pringle (1989), Armitage (2002) and Pessah, Chan 
& Psaltis (2008) have studied the role of magnetic fields within 
the boundary layer. 



the stellar magnctosphcre disrupts the inner accretion 
disk, leading to a magneto spheric mode of accretion in 
which gas becomes tied to stellar field lines and falls bal- 
listically on to the stellar surface (Konigl, 1991). The 
magnetic coupling between the star and its disk allows 
for angular momentum exchange, modifies the steady- 
state surface density profile close to the inner trunca- 
tion radius, and may allow the star to rotate more slowly 
than would otherwise be the case (Armitage & Clarke, 
1996; Collier Cameron & Campbell, 1993). Whether such 
"disk-locking" actually regulates the spin of young stars 
remains a matter of debate, however, and both theoret- 
ical and observational studies have returned somewhat 
ambiguous results (Herbst & Mundt, 2005; Matt & Pu- 
dritz, 2005; Rebull et al., 2006). 



3. Temperature profile 

Following the approach of Frank, King & Raine (2002), 
we derive the radial dependence of the effective temper- 
ature of an actively accreting disk by considering the net 
torque on a ring of width Ar. This torque - (dG/dr)Ar 
- does work at a rate, 



n—Ar = 
or 



d_ 

dr 



(GO) - GO' 



Ar 



(69) 



where O' = dO/dr. Written this way, we note that if we 
consider the whole disk (by integrating over r) the first 
term on the right-hand-side is determined solely by the 
boundary values of GO. We therefore identify this term 
with the transport of energy, associated with the viscous 
torque, through the annulus. The second term, on the 
other hand, represents the rate of loss of energy to the 
gas. We assume that this is ultimately converted into 
heat and radiated, so that the dissipation rate per unit 
surface area of the disk (allowing that the disk has two 
sides) is, 



D(r) 



GO' 9 
4?rr 8 



(70) 



where we have assumed a Keplerian angular velocity pro- 
file. For blackbody emission D(r) — o-T^ isk . Substituting 
for O, and for t/E using the steady-state solution given 
by equation (68), we obtain, 



3GM..M 



disk 




87rar 3 



(71) 



We note that, 



• Away from the boundaries (r >• r»), the tem- 
perature profile of an actively accreting disk is 
Idisk oc r~ 3 / 4 . This has the same form as for a 
passive disk given by equation (34). 

• The temperature profile does not depend upon the 
viscosity. This is an attractive feature of the the- 
ory given uncertainties regarding the origin and ef- 
ficiency of disk angular momentum transport. On 
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the flip side, it eliminates many possible routes to 
learning about the physics underlying v via obser- 
vations of steady-disks. 

Substituting a representative value for the accretion rate 



of M = 1(T 7 M fl yr" 



obtain for a Solar mass star at 



1 AU an effective temperature Xd; s k = 150 K. This is the 
surface temperature, as we will show shortly the central 
temperature is predicted to be substantially higher. 



4. Shakura-Sunyaev disks 

Molecular viscosity is negligible in protoplanetary 
disks. For a gas in which the mean free path is A, the 
viscosity 



Ac, 



(72) 



where c s is the sound speed. In turn, the mean free 
path is given by A = l/ner, where n is the number 
density of molecules with cross-section for collision a. 
These quantities are readily estimated. For example, 
consider a protoplanetary disk with £ = 10 3 gcm~ 2 and 
h/r = 0.05 at 1 AU. The midplane density is of the order 
of n ~ E/2mff/i ~4x 10 14 cm~ 3 , while the sound speed 
implied by the specified h/r is c s w 1.5 x 10 5 cms -1 . 
The collision cross-section of a hydrogen molecule is of 
the order of (Chapman & Cowling, 1970), 



u~2x 10~ 15 cm 2 , 



and hence we estimate, 



A <~ 1 cm 

~ 2x 10 5 cmV 1 . 



(73) 



(74) 



The implied disk evolution time scale r ~ r 2 /v then 
works out to be of the order of 10 13 yr - at least 10 6 
times too slow to account for observed disk evolution. 

In a classic paper, Shakura & Sunyaev (1973) noted 
that turbulence within the disk can provide an effec- 
tive viscosity that greatly exceeds molecular viscosity. 
For isotropic turbulence, the maximum scale of turbu- 
lent cells within the disk will be of the same order as the 
vertical scale height h, while the maximum velocity of 
turbulent motions relative to the mean flow is comparable 
to the sound speed c s (any larger velocity would lead to 
shocks and rapid dissipation of turbulent kinetic energy 
into heat). Motivated by such considerations, Shakura & 
Sunyaev (1973) proposed a parameterization, 



(75) 



where a is a dimensionless parameter that measures how 
efficient the turbulence is at creating angular momentum 
transport. We note at the outset that the existence of 
turbulence within the disk does not, a priori, guarantee 
that the outward angular momentum transport necessary 
to drive accretion will occur. 



In the standard theory of so-called "a-disks", a is 
treated as a constant. If this is done, it is possible to 
solve analytically for the approximate vertical structure 
of an actively accreting disk and derive a scaling for v as 
a function of r, E and a. Textbook discussions of this 
procedure can be found in Frank, King & Raine (2002), 
Armitage (2010), and many other places. Combining the 
known functional form for v with the disk evolution equa- 
tion (53) then yields a full theory for the predicted time 
dependence of the disk, with the only unknown being the 
appropriate value for a. This is all very well, but there is 
no physical reason to assume that a is a constant, and if 
instead a is regarded as a free function then much of the 
beguiling simplicity of the theory is lost, a-disk models 
should therefore be regarded as illustrative rather than 
definitive predictions for the evolution of the disk. 

It is straightforward to estimate how large a must be 
to account for the observed evolution of protoplanetary 
disks. Suppose, for example, that the evolution time scale 
at 50 AU is 1 Myr. Then starting from the a-prescription 
(equation 75), and noting that c s ~ htt, the evolution 
time scale becomes, 



r 
v 



1 



(76) 



Substituting for r and r, and assuming again that h/r = 
0.05, we obtain an estimate for a ~ 0.02. This is fairly 
typical - observational attempts to constrain a on large 
scales in protoplanetary disks (none of which are much 
more sophisticated than our crude estimate) tend to re- 
sult in estimates that are around 10~ 2 (Hartmann et al., 
1998) 6 . These values are an order of magnitude smaller 
than the values of a derived from the modeling of dwarf 
nova outbursts that occur in accretion disks around white 
dwarfs (Cannizzo, 1993; King, Pringle & Livio, 2007). 
It may be relevant to note that the disks around white 
dwarfs, and around other compact objects, are invariably 
more highly ionized than protoplanetary disks. 



5. Angular momentum transport processes 

Despite substantial recent progress, significant uncer- 
tainties persist as to the physical origin and properties 
of angular momentum transport within protoplanetary 
disks. The Reynolds number of the flow in the disk, 



Re ee 



UL 



(77) 



6 An important exception is modeling of the large-amplitude 
eruptive events known as FU Orionis outbursts (Hartmann & 
Kenyon, 1995), which, if interpreted as self-regulated thermal 
instabilities, require small values of a of the order of 10~ 3 or 
less (Bell & Lin, 1994). My own opinion is that these values are 
unreasonably small, and that FU Orionis events are instead ex- 
ternally triggered thermal instabilities that originate further out 
in the disk. 
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where U is a characteristic velocity and L a character- 
istic size scale, is extremely large (of the order of 10 14 
using the parameters that we previously estimated when 
considering the magnitude of molecular viscosity). Ter- 
restrial flows typically develop turbulence above a criti- 
cal Reynolds number of the order of 10 , so one's intu- 
ition would suggest that disk flows would surely be highly 
turbulent due to purely hydrodynamic effects. Detailed 
studies, however, do not support this conclusion. We first 
note that the condition for linear hydrodynamic stability 
in a differentially rotating fluid (the Rayleigh criterion) 
is that the specific angular momentum increase outward, 

(r 2 n) > 0. (78) 

dr 

In a Keplerian disk, r 2 f2 oc r 1 / 2 , so the flow is always 
linearly stable. 

A vast body of literature has investigated the possibil- 
ity of non-linear instabilities that might lead to turbu- 
lence within accretion disks. To date, there is no com- 
pelling evidence that such instabilities exist, and numer- 
ical simulations find that hydrodynamic perturbations 
in a Keplerian disk flow - which can in some circum- 
stances exhibit substantial transient growth (Afshordi, 
Mukhopadhyay & Narayan, 2005; Ioannou & Kakouris, 
2001) - ultimately decay (Balbus & Hawley, 2006; Bal- 
bus, Hawley & Stone, 1996; Shen, Stone & Gardiner, 
2006). Experiments yield a similar result (Ji et al., 2006). 
Although not all workers in the field would concur, the 
conclusion one draws from this is that in a model sys- 
tem - one in which the disk has simple thermodynamics 
(e.g. isothermal), is initially in an equilibrium state and 
is isolated from perturbing influences - hydrodynamic 
turbulence would not evolve spontaneously. 

Turbulence and angular momentum transport is pos- 
sible if one or both of magnetic fields and self-gravity are 
present. A sufficiently massive disk is unstable (Toomre, 
1964) to the development of trailing spiral arms, which 
act to transport angular momentum outward. We will 
discuss the physics underlying this instability later in the 
context of the disk instability model for giant planet for- 
mation, but for now we note that instability occurs when, 
roughly, 

^ > -• (79) 

Self-gravity could therefore play a role in protoplanetary 
disks at early epochs - when the disk may well be massive 
enough - but will not be important at late times when 
■^disk "C M*. Models for when self-gravity is important, 
and for the long-term evolution of disks evolving under 
the action of self-gravity, have been calculated by several 
authors (Clarke, 2009; Rafikov, 2009; Rice & Armitage, 
2009). The basic conclusion of such models is that - if 
other sources of angular momentum transport are weak 
or non-existent - then gas in the disk will settle into a 
stable self-gravitating state out to ~ 10 2 AU. Such disks 



are necessarily massive, and have a steep surface density 
profile. 

The hydrodynamic stability condition given by equa- 
tion (78) is also dramatically altered in the presence of 
even a weak magnetic field. Whereas a hydrodynamic 
flow is stable provided only that the specific angular 
momentum increase outward, a magnetohydrodynamic 
(MHD) flow requires that the angular velocity itself be 
an increasing function of radius, 

A (O 2 ) > 0, (80) 

in order to be stable (Balbus & Hawley, 1991; Chan- 
drasckhar, 1961; Vclikhov, 1959) 7 . This condition is 
not satisfied in Keplerian disks. As a consequence, in 
ideal (zero diffusivity) MHD an arbitrarily weak magnetic 
field suffices to render a Keplerian disk linearly unstable, 
with perturbations growing exponentially on a dynami- 
cal time scale. A comprehensive review of the physics of 
this instability - called the magnetorotational (MRI) or 
Balbus-Hawley instability - is given by Balbus & Haw- 
ley (1998). The MRI is a linear instability that leads to 
self-sustaining turbulence within sufficiently well-ionized 
accretion disks (Brandenburg et al., 1995; Stone et al., 
1996). It transports angular momentum outward, as is 
required to allow mass to flow inward and liberate gravi- 
tational potential energy. The magnitude of the effective 
a, generated by the MRI under ideal MHD conditions, 
has been estimated from local simulations to be of the 
order of a ~ 10" 2 (Davis, Stone & Pessah, 2010). This 
appears to be in encouraging agreement with the values 
inferred empirically for protoplanetary disks, but (as dis- 
cussed below) it must be remembered that ideal MHD is 
a poor approximation across much of the radial extent 
of real disks. Substantially more work is still required 
to determine the non-linear outcome of the MRI under 
realistic protoplanetary disk conditions. 

In addition to the relatively well-understood physics of 
the Rayleigh criterion, self-gravity, and the MRI, there 
are further possible sources of turbulence and angular 
momentum transport that have as their origin entropy 
gradients. The simplest to consider is convection in the 
vertical direction, which will occur in a disk (as in a star) 
if the vertical entropy profile is unstable. For many years 
it was thought that convection in disks transports angu- 
lar momentum inward, and that as a consequence it could 



7 The significance of Chandrasckhar's result for the origin of turbu- 
lence within the protoplanetary disk was appreciated by Safronov 
(1969), who noted that the MHD stability criterion does not re- 
duce to the Rayleigh criterion as the magnetic field tends toward 
zero, and that "for a weak magnetic field the cloud should be less 
stable than we found earlier in the absence of the field" . Safronov 
then, however, dismisses the MRI on the (incorrect) grounds that 
the instability requires that the viscosity and diffusivity are iden- 
tically zero. The importance of the MRI for accretion disks was 
only appreciated more than 20 years later by Balbus & Hawley. 
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not play any role whatsoever in disk evolution (Ryu & 
Goodman, 1992; Stone & Balbus, 1996). More recent 
simulations, however, have demonstrated that convec- 
tion does yield a positive value of a (Lcsur & Ogilvic, 
2010) 8 . Currently, it seems as if this result is primarily 
of academic interest, since it is unlikely that a sufficiently 
unstable vertical entropy profile can be sustained within 
a disk. Nonetheless, it is not impossible that there could 
be - at the least - zones within the disk where angular 
momentum transport via convection is non- negligible. 

In the radial direction, the condition for a rotating flow 
to be stable to linear axisymmetric disturbances in the 
presence of an entropy gradient is known as the Solberg- 
Hoi'land criterion. It can be written as, 

N 2 + n 2 >Q, (81) 

for a Keplerian disk, where N, the Brunt- Vaisala, fre- 
quency, is, 

7 p dr dr \pi ) ' v ' 

with 7 the adiabatic index. Protoplanetary disks are sta- 
ble to radial convection by this criterion. They can, how- 
ever, be unstable to a local, finite amplitude instability 
that is driven by the radial entropy gradient. This insta- 
bility, called the subcritical baroclinic instability (Lesur & 
Papaloizou, 2010; Petersen, Stewart & Julien, 2007), is 
present when, 

N 2 < 0, (83) 

(i.e. when the disk is Schwarzschild unstable), and there 
is either significant thermal diffusion or a thermal balance 
set by irradiation and radiative cooling. The subcritical 
baroclinic instability results in the formation of vortices 
(Klahr & Bodenhcimcr, 2003) - hydrodynamic structures 
of the type exemplified by Jupiter's Great Red Spot with 
non-zero vorticity to = V x v. Vortices are of particular 
interest because they can both transport angular momen- 
tum and, by trapping dust within their cores, accelerate 
the formation of larger solid bodies (Barge & Sommc- 
ria, 1995). How efficient they are at accomplishing these 
tasks is quite hard to assess, because in three dimensional 
disks vortices are subject to disruptive instabilities (Bar- 
ranco & Marcus, 2005; Lesur & Papaloizou, 2009; Lith- 
wick, 2009; Shen, Stone & Gardiner, 2006) 9 . At least 
qualitatively, the emerging picture is that the population 
of vortices present in a disk will reflect a balance between 



This flip is not due to a mistake on anyone's part. Rather, 
the early simulations appear to have had insufficient resolution, 
which led to a nearly axisymmetric pattern of convection that 
cannot transport angular momentum outward. 
In two dimensions, on the other hand, vortices are known to 
be long lived and quite effective agents of angular momentum 
transport (Godon & Livio, 1999; Johnson & Gammic, 2005). 



mechanisms that generate vorticity (such as the subcrit- 
ical baroclinic instability, which will probably work best 
when the disk is passive and not too optically thick) , and 
instabilities that destroy it. There is, as yet, no quanti- 
tative understanding of this balance, and hence it is not 
known whether vortices play a major role in either disk 
evolution or the early stages of planet formation. 

6. Layered disks 

The MRI is widely considered to be the most impor- 
tant (and possibly the only) source of angular momentum 
transport in accretion disks around white dwarfs, neutron 
stars and black holes. However in protoplanetary disks 
an interesting complication arises because the low ion- 
ization fraction leads to a finite conductivity. Resistivity 
(or other departures from ideal MHD due to ambipolar 
diffusion and the Hall effect) can then potentially damp 
the MRI and suppress turbulence and resulting angular 
momentum transport. The linear physics in this regime 
has been analyzed in numerous papers, including works 
by Blaes & Balbus (1994), Desch (2004) and Salmeron 
& Wardle (2005). A recent review by Balbus (2009) pro- 
vides a good entry into this literature. Here, I summarize 
a broad picture based on simple physics that was pro- 
posed by Gammie (1996) as a model for how the MRI 
might operate in protoplanetary disks. I stress that, in 
contrast to the prior discussion of the basic principles 
of the MRI, this application to the protoplanetary disk 
remains somewhat speculative. 

Following Gammie (1996), we begin by noting that in 
the presence of resistivity (assumed to be the most im- 
portant non-ideal MHD effect affecting the growth of the 
MRI) the magnetic field obeys the usual induction equa- 
tion, 

— = V x (v x B) - V x (f>V x B) , (84) 

where rj is the magnetic diffusivity. In turn, r\ can be 
written in terms of the electron fraction x = n e /nn via, 

n = 6.5 x lO" 3 ^ 1 cmV 1 . (85) 

Our goal is to determine the minimum x for which the 
MRI will be able to operate despite the damping caused 
by the diffusivity. To do this, we note that resistivity 
damps small scales most easily. We therefore consider 
the largest disk scale I = h, and equate the MRI growth 
time scale (Balbus & Hawley, 1998), 

h , s 

TMRI ~ (86) 

VA 

where va = \J B 2 / {An p) is the Alfven speed, with the 
damping time scale, 

h 2 

Tdamp ~ ■ (87) 
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FIG. 15 Schematic illustration of the layered disk model for 
protoplanetary disks. In this model, the innermost regions 
of the disk are hot enough that thermal ionization suffices to 
couple the magnetic field to the gas and allow the MRI to 
operate. At large radii, some combination of stellar X-rays 
and cosmic rays penetrate the entire thickness of the disk and 
provide the necessary ionization. At intermediate radii, it is 
hypothesized that accretion occurs primarily in an active sur- 
face layer of whose column depth is set by the penetration 
power of the dominant non-thermal ionization process (either 
cosmic rays, as in the original version of this model proposed 
by Gammie (1996), or stellar X-rays in some more recent ver- 
sions). Underlying this layer is a magnetically inactive "dead 
zone" in which the MRI is suppressed and the level of turbu- 
lence is greatly weakened. 

This yields a simple criterion for the MRI to operate: 

T] < hv A . (88) 

It remains to estimate appropriate values for h and va- 
For a crude estimate, we can guess that at 1 AU h ~ 
10 12 cm and that va ~ c s ~ 10 5 cms -1 (more accurately, 
va ~ a 1 / 2 c s in MRI turbulence that yields an effective 
Shakura-Sunyaev a). In that case the limit becomes r\ < 
10 17 cm 2 s _1 which translates into a minimum electron 
fraction, 

x > 1CT 13 , (89) 

which is more or less the 'right' value derived from 
more rigorous analyses (Balbus & Hawley, 1998; Gam- 
mie, 1996). The most important thing to note is that 
this is an extremely small electron fraction! The linear 
MRI growth rate is so large that a tiny electron fraction 
couples the gas to the magnetic field well enough that the 
MRI can overcome the stabilizing influence of diffusion. 

Although only a small degree of ionization is required 
for the MRI to work, there may be regions in the pro- 
toplanetary disk where even x ~ 10~ 13 is not attained. 
Considering first thermal ionization processes, calcula- 
tions of collisional ionization by Umebayashi (1983) show 
that ionization of the alkali metals suffices to drive x > 
1CP 13 . This, however, requires temperatures T ps 10 3 K 
and above. Only the very innermost disk - within a few 
tenths of an AU of the star - will therefore be able to 
sustain the MRI as a result of purely thermal ionization. 

At larger disk radii the ionization fraction will be de- 
termined by a balance between non-thermal ionization 



processes and recombination. Two sources of ionization 
are of particular interest, 

• Ionization by stellar X-rays. T Tauri stars are ob- 
served to be strong X-ray sources (Feigelson et al., 
2007) , and the harder components of the emission 
are unquestionably penetrating enough to ionize a 
fraction of the column through the disk. 

• Ionization by cosmic rays. Cosmic rays have a 
stopping length that is of the order of Si ayC r = 
100 gem" 2 (Umebayashi & Nakano, 1981). If 
present they are therefore likely to be more pen- 
etrating and important than X-rays. It is possible 
(or perhaps probable), however, that the disk may 
be screened from the interstellar cosmic ray flux 
by the magnetized plasma flowing away from the 
system in a wind. 

If no other sources of ionization exist, then at radii where 
the disk is simultaneously too cool to be collisionally ion- 
ized, and dense enough that the interior is shielded from 
X-rays (or cosmic rays, if present) we might expect a 
novel structure in which the disk near the midplane is 
quiescent (a "dead zone") and only a surface layer is 
magnetically active and supporting accretion (Gammie, 
1996). Such a layered disk model is depicted in Figure 15. 

Layered disk models are qualitatively distinct from or- 
dinary fully viscous disks in two critical respects. First, 
the mass flux through the active layer is set (by analogy 
to equation 68) by the product of the viscosity and col- 
umn z4 a yerSi a yer in the active layer, which is independent 
of the total column density at that radius. The amount 
of mass that the active layer can support is predicted to 
be a decreasing function of radius, so gas flowing inward 
from large radii "drops out" of the flow and accumulates 
in the dead zone. Second, layered models predict that the 
midplane of the disk ought to be almost quiescent at pre- 
cisely those radii of greatest interest for planet formation. 
This has important implications for the settling of dust, 
for the subsequent growth of planetesimals, and for the 
migration of low mass planets (Matsumura & Pudritz, 
2005). There are also possible implications for variabil- 
ity, since the growing dead zone provides a reservoir of 
gas at small radii which could in principle be heated and 
activated leading to a burst of accretion (Armitage, Livio 
& Pringle, 2001; Zhu, Hartmann & Gammie, 2009). 

Given these intriguing possibilities, do layered disks ac- 
tually exist? Currently there are few observational con- 
straints, though theoretically the idea remains plausible. 
Numerical simulations (Sano & Stone, 2002) have con- 
firmed that when the magnetic Reynolds number (eval- 
uated on the relevant scale for a disk), 

Re M = V \ (90) 

falls below a critical value that is in the range of 1 
30, turbulence driven by the MRI is suppressed. This is 
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broadly consistent (i.e. it yields a similar critical elec- 
tron fraction) with the simple arguments outlined above. 
However, very little energy would be required to maintain 
a high enough level of ionization for the MRI to operate, 
and it has been suggested that enough power from the 
turbulence could couple into the ionization to sustain an 
MRI-active disk (Inutsuka & Sano, 2005). Another pos- 
sibility — which receives support from recent numerical 
simulations by Turner, Sano & Dziourkevitch (2007) — is 
that turbulent mixing might be able to transport enough 
charge from the surface layers to the midplane to allow 
non-zero transport to occur there. Given these theoreti- 
cal uncertainties, the question of whether protoplanetary 
disks have a structure akin to that shown in Figure 15 
remains open. A critical question is whether the rate of 
recombination within the gas is high enough to mop up 
free electrons efficiently Recombination depends upon 
difficult aspects of the chemistry and dust physics within 
the disk (e.g. how many small dust particles are present 
close to the disk surface, what is the abundance of metal 
ions in the gas phase?), and the resulting uncertainties 
are at least as large as those arising from the nature of 
the ionization or from the basic physics of the MRI under 
non-ideal conditions. 



7. Disk dispersal 

Loss of the gaseous component of protoplanetary disks 
sets a time limit for the completion of gas giant forma- 
tion, and will affect the environment for terrestrial planet 
formation as well. The self-similar solution for the evo- 
lution of a viscous disk (equation 62) predicts that the 
surface density and accretion rate decline as power-laws 
at late times, and hence that the transition between disk 
and diskless states should be gradual. Observationally, 
this is not what is observed. Young stars with proper- 
ties intermediate between CTTS and WTTS (so called 
"transition objects", though this class has multiple def- 
initions in the literature) are relatively uncommon, and 
make up of the order of 10% of the total population of 
stars that display at least one signature of a circumstel- 
lar disk. From these statistics, one infers that the time 
scale for clearing the disk is short - of the order of 10 5 yr 
(Simon & Prato, 1995; Wolk & Walter, 1996). 

A clue to the mechanism that may drive disk disper- 
sal was provided by HST observations of low mass stars 
exposed to the strong ionizing flux produced by massive 
stars in the core of the Orion Nebula's Trapezium cluster 
(O'Dell, Wen & Hu, 1993). The images reveal tadpole- 
shaped nebulae surrounding young stars with circumstcl- 
lar disks, which are interpreted as the signature of pho- 
toevaporation and escape of disk gas as a result of illumi- 
nation by external ionizing radiation (Johnstone, Hollen- 
bach & Bally, 1998). The physics underlying this process 
is relatively simple (Bally & Scoville, 1982; Hollcnbach 
et al., 1994; Shu, Johnstone & Hollenbach, 1993), and 
is closely related to the well-studied problem of Comp- 



ton heated winds from accretion disks in Active Galactic 
Nuclei (Begelman, McKee & Shields, 1983). Extreme ul- 
traviolet (EUV) photons with E > 13.6 eV ionize and 
heat a surface layer of the disk, raising it to a temper- 
ature T ~ 10 4 K characteristic of an HII region. The 
sound speed in the photoionized gas is c s ~ 10 kms -1 . 
Outside a critical radius r g , given by, 



the sound speed in the hot gas exceeds the local Keplerian 
speed. The gas is then unbound, and flows away from 
the disk as a thermal wind. For a Solar mass star, r g as 
estimated by equation (91) is at 9 AU. 

The same basic process can occur regardless of whether 
the EUV flux arises from an external source, such as a 
massive star in a cluster, or from the central star itself. 
In the typical star formation environment (Lada & Lada, 
2003), however, most low mass stars receive too low a 
dose of EUV radiation from external sources to destroy 
their disks (Adams et al., 2006). Photoevaporation due 
to radiation from the central star is therefore likely to be 
necessary for disk dispersal. In this regime, Hollcnbach 
et al. (1994) derived an estimate for the mass loss rate 
due to photoevaporation, 

4 * 10 - w {wh*)" 2 {^)" 2 m^'- 1 

(92) 

where $ is the stellar ionizing flux. Most of the wind mass 
loss is predicted to originate close to r g , with a radial de- 
pendence of the mass loss given by £ cx r~ 5 / 2 . Numerical 
hydrodynamic simulations by Font et al. (2004) largely 
confirm this basic picture, although in detail it is found 
both that mass is lost for radii r < r g and that the in- 
tegrated mass loss is a factor of a few smaller than that 
predicted by the above equation. 

The combination of a photoevaporative wind and vis- 
cous disk evolution leads to rapid disk dispersal (Clarke, 
Gendrin & Sotomayor, 2001). Calculations by Alexan- 
der, Clarke & Pringle (2006) suggest a three-stage sce- 
nario depicted schematically in Figure 16, 

• Initially M ^> M win( j. The wind mass loss has neg- 
ligible effect on the disk, which evolves in a similar 
way to an ordinary viscous model. The mass accre- 
tion rate and surface density gradually drop on the 
viscous time scale of the entire disk (determined at 
large radii), which is of the order of a Myr. 

• After a few Myr, the accretion rate drops suffi- 
ciently so that M ~ Mwi n d. The wind is then 
strong enough to dominate the disk surface den- 
sity evolution near r g , opening up a gap in the disk 
and cutting off the inner disk from resupply by gas 
flowing in from the reservoir at larger radii. The 
inner disk then drains on to the star on its own 
(short) viscous time scale, which can be of the or- 
der of 10 5 yr or less. 
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FIG. 16 Schematic depiction of how photoevaporation driven 
by a central source of UV radiation is predicted to disperse 
the protoplanetary disk. In the initial phase, UV radiation 
(shown as the red arrows) ionizes the surface of the disk, pro- 
ducing a vertically extended bound atmosphere for r < r g and 
mass loss in a thermal wind for r > r g . The ionizing flux that 
photoevaporates the outer disk arises primarily from stellar 
photons scattered by the atmosphere at small radii (the 'dif- 
fuse field'). After several Myr, the disk accretion rate drops 
to a value that is of the same order as the wind mass loss 
rate. At this point, the wind opens up a gap in the disk close 
to r g , cutting off the inner disk from resupply by the disk 
further out. The inner disk then drains rapidly on to the star 
- producing an inner hole - and the direct UV flux from the 
star photoevaporates the outer region. 



• Once the inner disk has drained, the remaining gas 
in the outer disk is directly illuminated by UV radi- 
ation from the star (previously, the dominant flux 
was photons scattered on to the outer disk from a 
bound atmosphere at smaller radii) . This radiation 
rapidly burns through the outer disk removing all 
remaining gas. 

The primary source of uncertainty in these models is the 
origin and magnitude of the stellar ionizing flux. There 
are few constraints on $ for Solar mass T Tauri stars 
(Alexander, Clarke & Pringle, 2005), and essentially no 
information on any scaling with stellar mass. 

To date EUV-driven photoevaporation models for disk 
evolution have received the most theoretical attention. 
This, however, is primarily because the physics of EUV- 
ionized gas is particularly easy to calculate. Qualita- 
tively similar flows can be driven by softer FUV radia- 



tion (6 eV < E < 13.6 eV), which suffices to dissociate 
H2 molecules and drives evaporative flow from the outer 
disk where the escape velocity is smaller. The detailed 
physics of such flows - which resemble photodissociation 
regions rather than HII regions - is harder to calculate 
because the temperature of the heated gas is determined 
by a balance between grain photoelectric heating and 
cooling by both atomic and molecular lines (Adams et 
al., 2006; Gorti & Hollenbach, 2009). Harder X-ray pho- 
tons can also be important, with recent work suggesting 
that X-rays may in fact dominate the total photoevapo- 
rative mass loss rate for protoplanetary disks (Ercolano, 
Clarke & Drake, 2009). 



D. The condensation sequence 

In an actively accreting disk, there must be a temper- 
ature gradient in z in order for energy to be transported 
from the dense midplane where it is probably liberated 
to the photosphere where it is radiated (note that for a 
thin disk with fc/r « 1 gradients in z will dominate over 
radial gradients, which can consistently be ignored). A 
simple application of the theory of radiative transport 
in plane-parallel media (Rybicki & Lightman, 1979) al- 
lows us to derive the relation between the central disk 
temperature T c and the effective disk temperature T^isk- 

To proceed, we define the optical depth to the disk 
midplane, 



(93) 



where kr is the Rosseland mean opacity and £ is the 
disk surface density. The vertical density profile of the 
disk is p(z). If the vertical energy transport occurs via 
radiative diffusion (in some regions convection may also 
be important), then for r 3> 1 the vertical energy flux 
F{z) is given by the equation of radiative diffusion, 



F(z) = 



WaT 3 dT 
3krp dz 



(94) 



Let us assume for simplicity that all the energy dissipa- 
tion occurs at z = 10 . In that case F(z) — cT^isk is 
a constant with height. Integrating assuming that the 
opacity is a constant, 
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(95) 



10 Note that for disks in which the MRI is active numerical simula- 
tions by Miller & Stone (2000) suggest that a significant fraction 
of the liberated energy is transported to high z and dissipated at 
much smaller optical depths, possibly forming a hot "corona". 
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where for the final equality we have used the fact that 
for t ^s> 1 almost all of the disk gas lies below the photo- 
sphere. For large r we expect that T c 4 » T^ isk , and the 
equation simplifies to, 



TABLE II Condensation temperatures for selected materials 
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J disk 



(96) 



The implication of this result is that active disks with 
large optical depths are substantially hotter at the mid- 
plane than at the surface. For example, if r = 10 2 to 
the thermal radiation emitted by the disk at some radius 
then T c sa 3Tdi s k- This is important since it is the central 
temperature that determines the sound speed that enters 
into the viscosity (equation 75), and it is also the central 
temperature that determines which ices or minerals can 
be present. Relatively modest levels of accretion can thus 
affect the thermal structure of the disk substantially. 

For both terrestrial planet formation, and gas giant 
planet formation if it occurs via the core accretion mech- 
anism, the evolution of the trace solid component of the 
disk is of great interest. The gas that forms the proto- 
planetary disk will contain interstellar dust grains made 
up of a mixture of silicates, graphite and polycyclic aro- 
matic hydrocarbons (PAHs) . In the interstellar medium, 
measurements of the wavelength dependence of extinc- 
tion can be fit by assuming that the dust grains follow a 
power-law size distribution (Mathis, Rumpl & Nordsicck, 
1977), 



n(a) oc a 



-3.5 



(97) 



where a is the grain size (assumed to be spherical) and 
the distribution extends from 0.005 /im to about 1 /ini. 
This distribution is generally assumed to be the starting 
point for further evolution within the denser conditions 
prevailing within the disk. In the hottest, inner regions 
of the disk the central temperature can be high enough 
to destroy the grains (1000 K to 2000 K, depending on 
whether the grains are made of carbon or silicate). The 
resulting absence of dust very close to the star consti- 
tutes one of the main arguments against an in situ origin 
for hot Jupiters, but the dust destruction radius is suffi- 
ciently close in (normally substantially less than 1 AU) 
that it rarely impacts either terrestrial or, especially, gi- 
ant planet formation. It is, however, important observa- 
tionally, since once the dust is destroyed the remaining 
gas phase opacity is greatly reduced. There will there- 
fore be an opacity "hole" in the inner disk even if gas is 
present there. 

If the gas that makes up the protoplanetary disk has 
a known elemental composition (for example that of the 
Sun), then it is a well defined problem (for a chemist!) 
to calculate the most thermodynamically stable mix of 
chemical species at any given pressure and temperature. 
The abundance of different minerals and ices within the 
disk will follow this condensation sequence provided that 
there is sufficient time for chemical reactions to reach 
equilibrium - this may be a reasonable assumption in 
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the hot inner disk but deviations will occur due to slow 
chemical reactions in the cool outer disk and radial drift 
of both gas and solids. The equilibrium mix depends 
more strongly on temperature than on pressure, so we 
can roughly map the condensation sequence into a pre- 
dicted variation of disk composition with radius. Ta- 
ble II, adapted from Lodders (2003), lists characteris- 
tic temperatures below which different ices and min- 
erals are predicted to be dominant. Of these, by far 
the most important is the temperature below which wa- 
ter ice can be present - this is 180 K at a pressure of 
10~ 4 bar (though for the conditions in the protoplane- 
tary disk, water ice requires moderately cooler conditions 
with T « 150 K). For a Solar mix of elements, the sur- 
face density of condensible materials rises dramatically 
once water ice forms, 



S(ices + rock) ~ 4E(rock) 



(98) 



though the ratio depends upon still uncertain determi- 
nations of the exact Solar composition (Lodders, 2003). 
It is tempting - and extremely plausible - to associate 
changes in the efficiency or outcome of planet formation 
(in particular the division between terrestrial and gas gi- 
ant planets in the Solar System) with the large change in 
the predicted surface density of solids that occurs at this 
radius. 

The radius in the protoplanetary disk beyond which 
water ice can be present is called the snow line. In the 
Solar System, water-rich asteroids are found only in the 
outer asteroid belt (Morbidelli et al., 2000), which sug- 
gests that the snow line in the Solar Nebula fell at around 
3 AU. Passive protoplanetary disks are predicted to have 
snow lines at substantially smaller radii - in some cases 
interior to 1 AU - though accretion rates within the plau- 
sible range for Classical T Tauri disks suffice to push the 
snow line out to the inferred location of 3 AU (Garaud 
& Lin, 2007; Lecar et al., 2006). 



III. PLANET FORMATION 

The formation of planets from sub-micron size dust 
particles requires growth through at least 12 orders of 
magnitude in spatial scale. It is useful to consider dif- 
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ferent size regimes in which the interaction between the 
solid component and the gas is qualitatively distinct: 

• Dust - small particles ranging from sub-micron to 
cm in scale. These particles are well-coupled to the 
gas, but there can be slow drift either vertically or 
radially. Growth occurs through physical collisions 
leading to agglomeration. 

• "Rocks" - objects of meter scale. These particles 
have increasingly weak coupling to the gas, and so 
it can be useful to approximate their dynamics as 
being a combination of Keplerian orbits plus aero- 
dynamic drag. Growth through this size regime is 
deduced to be rapid but the mechanism remains 
uncertain. 

• Planetesimals - bodies of 10 km scale and above. 
Planctesimals are massive enough that their dy- 
namics is largely decoupled from that of the gas. 
A population of bodies of this size is often consid- 
ered as the initial condition for subsequent stage of 
planet formation, since the evolution of such a pop- 
ulation is a well-posed - albeit difficult - N-body 
problem involving primarily purely gravitational 
forces (though for smaller planetesimals, questions 
regarding the bodies material strength can also be 
pertinent). 

• Earth mass planets or progenitors of the giant 
planet cores. Once growing planets reach masses 
of the order of that of Earth, they again become 
coupled to the gas disk, though this time via gravi- 
tational rather than aerodynamic interactions. We 
will discuss this coupling later in the context of dif- 
ferent regimes of planetary migration. For terres- 
trial planet formation it is possible that the forma- 
tion of Earth mass bodies occurs after the gas disk 
has been dispersed (in which case this coupling is 
moot), but for growing giant planet cores it is in- 
evitable that interaction will take place. 

• Planetary cores with masses of the order of 
10 M®. At around this mass, there is a transition 
from a quasi-hydrostatic core + envelope structure 
to a regime of rapid gas accretion. 

Although it predates the discovery of extrasolar plane- 
tary systems, the review by Lissauer (1993) remains an 
excellent, readable summary of much of the physics that 
we will address in this section. 



A. Planetesimal formation 

A spherical particle of radius a, moving relative to the 
gas at velocity v, experiences an aerodynamic drag force 
Fd that opposes its motion, 



where Cd is the drag coefficient. The form of the drag co- 
efficient depends upon the size of the particle compared 
to the mean free path A of molecules in the gas (Wei- 
denschilling, 1977b; Whipple, 1972). For small particles 
(typically of cm size or less) for which, 



a < -A 

4 



the drag coefficient is, 



C D = 



3v 



(100) 



(101) 



where v = (8/7r) 1 / 2 c s is the mean thermal velocity in 
the gas. This is called the Epstein regime of drag. For 
larger particles the Stokes drag law is valid. Defining the 
Reynolds number via, 



Re = 



2av 



(102) 



where v is the microscopic (molecular) viscosity in the 
gas, the drag coefficient can be expressed as a piecewise 
function, 



C D = 24RC" 1 Re < 1 

C D = 24Re-°' 6 1 < Re < 800 

C D = 0.44 Re > 800. 



(103) 



We will apply these drag laws to consider both the ver- 
tical distribution and radial drift of small solid bodies 
within the gas disk. 



1. Dust settling 

Dust particles are strongly coupled to the gas via drag 
forces. For a particle of mass m, the friction time scale 
is defined as, 



tfr 



Fd 



(104) 



It is the time scale on which drag will lead to order unity 
changes in the relative velocity between the particle and 
the gas. Writing the particle mass m = (4/3)-7ra 3 p (i in 
terms of the material density pa, the friction time scale 
has a simple form in the Epstein regime, 



f-tric — _ • 
p V 



(105) 



Adopting conditions appropriate to 1 AU within the disk, 
p = 5 x 10~ 10 gem -3 , v = 2.4 x 10 5 cms -1 and p d = 
3 gem -3 we obtain if r ; c w 2.5 s. Small particles are thus 
very tightly coupled to the gas. 

Consider a thin, vertically isothermal gas disk with sur- 
face density £ and scale height h = c s /Qk- The vertical 
density profile is, 



F D = ~\c D ■ na 2 ■ pv 2 



(99) 



p{z) = 



-z 2 /2h 2 



(106) 
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To start with, let us ignore the effects of turbulence and 
assume that the disk is entirely quiescent. In this case the 
important forces acting on a particle at height z above 
the midplane are the vertical component of gravity and 
gas drag, given by, 
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Given the strong coupling expected for dust particles ter- 
minal velocity will rapidly be attained, so we equate these 
to obtain the settling speed, 
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(108) 



Settling is more rapid at higher z (where the gas density 
is lower and the vertical component of gravity stronger) , 
and for larger grains. For example, for micron sized dust 
particles at z = h at 1 AU the settling velocity is w se ttie ~ 
0.1 cms -1 and the settling time scale, 



^settle 



f settle 



- 2 x 10 5 yr. 



(109) 



In the absence of turbulence, then, we expect micron 
sized dust particles to sediment out of the upper layers 
of the disk on a time scale that is short compared to the 
disk lifetime, while for particles with sizes < 0.1 pm the 
time scale is marginal. 

Only the density in the equation for the settling time 
scale is a function of height. Inserting the expression for 
the vertical density profile the general expression for the 
settling time scale becomes, 



(2006), present detailed calculations of the impact of tur- 
bulence (primarily that generated by the MM) on parti- 
cle settling. For a simple treatment, that is at least qual- 
itatively correct, we can assume that the effective dust 
diffusion coefficient in z is simply equal to the turbulent 
viscosity, 



v<i = v = ac s h 



(111) 



and that the time scale over which turbulence will stir 
up a layer of thickness z is given by the usual diffusive 
expression, 



^stir — 

Vd 



(112) 



Clearly i st ; r increases with z, whereas i se ttie decreases 
with height. Equating these expressions, we obtain an 
estimate of the thickness of the dust layer that will re- 
main well-mixed as a consequence of turbulence. The 
thickness z is given implicitly via, 



-z 2 /h 2 _ ^PdO. 

2aS ' 



(113) 



Substituting typical disk and dust properties on the right 
hand side of this expression, we find that for micron sized 
particles in a disk with a = 10~ 2 and S = 10 3 gem -2 
the right hand side is of the order of 10 -4 . This means 
that small dust particles will remain suspended in the 
disk up to at least several scale heights. On the other 
hand, for a — 1 cm the right hand side is around unity 
- so large particles will settle to z ~ h even if the disk is 
fully turbulent. 
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(110) 2. Settling with coagulation 



The strong z-dependence implies that dust will settle out 
of the upper regions of the disk rather rapidly in the ab- 
sence of turbulence. This is of some interest since scat- 
tered light images of protoplanetary disks (e.g. Burrows 
et al., 1996) are sensitive to dust well away from the mid- 
plane. 

Assuming that the disk is quiescent is not very realistic 
- the same turbulence that leads to angular momentum 
transport is very likely to "stir up" the dust and prevent 
it settling to the midplane as easily as our estimate sug- 
gests. Working out the influence of turbulence on particle 
settling is quite tricky, since it involves two subtle issues, 

• What is the effective diffusion coefficient in the ver- 
tical direction for passive tracer particles injected 
into the turbulent flow? 

• How well coupled are the particles to the gas? 

Several authors, including Dullemond & Dominik (2004), 
Johanscn & Klahr (2005), Carballido, Stone & Pringle 
(2005), Turner et al. (2006) and Fromang & Papaloizou 



Even if the neglect of turbulence was justified - and 
it is not - the estimate of the dust settling time in a 
laminar would be incomplete because it ignores the like- 
lihood that dust particles will collide with one another 
and grow during the settling process. The settling ve- 
locity increases with the particle size, so any such coag- 
ulation hastens the collapse of the dust toward the disk 
midplane. 

To estimate how fast particles could grow during sed- 
imentation we appeal to a simple single particle growth 
model (Dullemond & Dominik, 2005; Safronov, 1969). 
Imagine that a single "large" particle, of radius a and 
mass to = (4/3)na 3 pd, is settling toward the disk mid- 
plane at velocity u se ttie through a background of much 
smaller solid particles. By virtue of their small size, the 
settling of the small particles can be neglected. If every 
collision leads to coagulation, the large particle grows in 
mass at a rate that reflects the amount of solid material 
in the volume swept out by its geometric cross-section, 
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FIG. 17 The settling and growth of a single particle in a lami- 
nar (non-turbulent) protoplanetary disk. The model assumes 
that a single particle (with initial size a = 1 pm (solid line), 
0.1 Jim (dashed line), or 0.01 pm (long dashed line) accretes 
all smaller particles it encounters as it settles toward the disk 
midplane. The smaller particles are assumed to be at rest. 
The upper panel shows the height above the midplane as a 
function of time, the lower panel the particle radius a. For 
this example the disk parameters adopted are: orbital radius 
r = 1 AU, scale height h = 3 x 10 cm, surface density 
E = 10 3 g cm" 2 , dust to gas ratio / = 10~ 2 , and mean ther- 
mal speed v = 10 5 cm s _1 . The dust particle is taken to have 
a material density pd = 3 g cm~ 3 and to start settling from 
a height zo = 5h. 

where / is the dust to gas ratio in the disk. Substituting 
for the settling velocity one finds, 

dm Ztff 

— — = zm. 115 

At 4 v v ' 

Since z = z(t) this Equation cannot generally be inte- 
grated immediately 11 , but rather must be solved in con- 
cert with the equation for the height of the particle above 
the midplane, 

^=-^V, (116) 
At p v v ' 

Solutions to these equations provide a very simple model 
for particle growth and sedimentation in a non-turbulent 
disk. 



Note however that if the particle grows rapidly (i.e. more rapidly 
than it sediments) then the form of the equation implies expo- 
nential growth of m with time. 



Figure 17 shows solutions to equations (115) and (116) 
for initial particle sizes of 0.01 /zm, 0.1 /zm and 1 pm. The 
particles settle from an initial height zq = 5h through a 
disk whose parameters are chosen to be roughly appro- 
priate to a (laminar) Solar Nebula model at 1 AU from 
the Sun. Both particle and growth and vertical settling 
are extremely rapid. With the inclusion of coagulation, 
particles settle to the disk midplane on a time scale of 
the order of 10 3 yr - more than two orders of magnitude 
faster than the equivalent time scale in the absence of 
particle growth. By the time that the particles reach the 
midplane they have grown to a final size of a few mm, 
irrespective of their initial radius. 

The single particle model described above is very sim- 
ple, both in its neglect of turbulence and because it as- 
sumes that the only reason that particle-particle colli- 
sions occur is because the particles have different ver- 
tical settling velocities. Other drivers of collisions in- 
clude Brownian motion, turbulence, and differential ra- 
dial velocities. The basic result, however, is confirmed 
by more sophisticated models (Dullemond & Dominik, 
2005), which show that if collisions lead to particle adhe- 
sion growth from sub-micron scales up to small macro- 
scopic scales (of the order of a mm) occurs rapidly. There 
are no time scale problems involved with the very earliest 
phases of particle growth. Indeed, what is more problem- 
atic is to understand how the population of small grains 
- which are unquestionably present given the IR excesses 
characteristic of Classical T Tauri star - survive to late 
times. The likely solution to this quandary involves the 
inclusion of particle fragmentation in sufficiently ener- 
getic collisions, which allows a broad distribution of par- 
ticle sizes to survive out to late times. Fragmentation 
is not likely given collisions at relative velocities of the 
order of a cm s _1 - values typical of settling for micron- 
sized particles - but becomes more probable for collisions 
at velocities of a m s _1 or higher. 

3. Radial drift of particles 

Previously we showed (equation 48) that the azimuthal 
velocity of gas within a geometrically thin disk is close 
to the Keplerian velocity. That it is not identical, how- 
ever, turns out to have important consequences for the 
evolution of small solid bodies within the disk (Weiden- 
schilling, 1977b). We can distinguish two regimes, 

• Small particles (a < cm) are well-coupled to the 
gas. To a first approximation we can imagine that 
they orbit with the gas velocity. Since they don't 
experience the same radial pressure gradient as the 
gas, however, this means that they feel a net in- 
ward force and drift inward at their radial terminal 
velocity. 

• Rocks (a > m) are less strongly coupled to the gas. 
To a first approximation we can imagine that they 
orbit with the Keplerian velocity. This is faster 
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than the gas velocity, so the rocks see a headwind 
that saps their angular momentum and causes them 
to spiral in toward the star. 

To quantify these effects, we first compute the magnitude 
of the deviation between the gas and Keplerian orbital 
velocities. Starting from the radial component of the 
momentum equation, 



GM, 



ldP 

p dr ' 



(117) 



we write the variation of the midplane pressure with ra- 
dius as a power-law near radius r , 



P = Pn 



ro 



where Pq = poc 2 s . Substituting, we find, 



«0,gas = V K (1 - V) 



1/2 



where 



V 



n- 
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(119) 



(120) 



'K 



Typically n is positive (i.e. the pressure decreases out- 
ward), so the gas orbits slightly slower than the local 
Keplerian velocity. For example, for a disk of constant 
h(r)/r = 0.05 and surface density profile S oc r _1 we 
have n = 3 and, 



v 4>,i 



0.996wx- 



(121) 



The fractional difference between the gas and Keplerian 
velocities is small indeed! However, at 1 AU even this 
small fractional difference amounts to a relative velocity 
of the order of 100 ms _1 . Large rocks will then experience 
a substantial, albeit subsonic, headwind. 

The effect of the drag force on the dynamics of particles 
of arbitrary sizes has been calculated by Wcidcnschilling 
(1977b). Here, we adopt the approach of Takeuchi & 
Lin (2002) and proceed by considering the radial and 
azimuthal equations of motion for the particle 12 , 
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We simplify the azimuthal equation by noting that the 
specific angular momentum always remains close to Ke- 
plerian (i.e. the particle spirals in through a succession 
of almost circular, almost Keplerian orbits), 

~ iv ^ (rv K ) = ^v r v K . (123) 



12 Although this calculation is straightforward, it's easy to confuse 
the three different azimuthal velocities that are involved — that 
of the particle, that of the gas, and the Kepler speed. Be careful! 



This yields, 
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Turning now to the radial equation, we substitute for 
using equation (119). Retaining only the lowest order 
terms, 
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The dv r /dt term is negligible, and for simplicity we also 
assume that v r ,gas <C v r , which will be true for those par- 
ticles experiencing the most rapid orbital decay. Elimi- 
nating (v^ — v<p,gas) between equations (124) and (125) 
we obtain, 
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This result can be cast into a more intuitive form by 
defining a dimensionless stopping time, 

Tfric = t{ T i c Q.K , (127) 

in terms of which the particle radial velocity is, 



v r 
v K 
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The peak radial velocity is attained when Tf r j c = 1 (i.e. 
when the friction time scale equals , and equals 
independent of the disk properties. 

Figure 18 plots v t /vk as a function of the dimension- 
less stopping time for a fiducial disk with h/r = 0.05. 
Using equations (101) and (103), one can associate a par- 
ticular Tf r i c with a unique particle size o given known con- 
ditions in the protoplanetary disk. Generically, one finds 
that at radii of a few AU the peak inspiral rate is attained 
for particles with size of the order of 10 cm to a few m. 
The minimum inspiral time scale at a given orbital radius 
depends only on 77 - at 1 AU it is of the order of 100 yr. 
The inescapable conclusion is that the radial drift time 
scale <C disk lifetime for meter-scale bodies in the 
protoplanetary disk. 

As we noted earlier, the fact that most of the heavy 
elements in the Solar System are found in the Sun means 
that we can tolerate some loss of planetary raw mate- 
rial during planet formation. However, radial drift time 
scales as short as 100 yr would clearly lead to a catas- 
trophic loss of mass into the star unless, in fact, growth 
through the meter-scale size regime is very fast. The 
most important conclusion from this analysis is, there- 
fore, that planetesimal formation must be a rapid pro- 
cess. This is a robust inference since it derives directly 
from the unavoidable existence of a velocity differential 
between the gas disk and solid bodies orbiting within it. 

The radial drift velocities given by equation (128) im- 
ply significant radial migration over the lifetime of the 
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FIG. 18 Radial drift velocity of particles at the midplane of 
a protoplanetary disk with h/r = 0.05, plotted as a function 
of the dimensionless stopping time Tf r i c . The radial velocity 
of the gas has been set to zero. The most rapid inward drift 
occurs for a physical stopping time Q.^ 1 , which for typical disk 
models translates to a particle size in the 10 cm to m range. 
At 1 AU, the peak inward velocity is around 60 ms _1 , which 
implies a decay time of less than 100 yr. 



disk - not just for particles at the most vulnerable meter- 
scale size range but also for substantially smaller and 
larger bodies. This means that we should expect substan- 
tial changes in the local ratio of solids to gas as a function 
of time and radius in the disk (Takeuchi, Clarke & Lin, 
2005). Under some circumstances, radial drift may allow 
solids to pileup within the inner disk, potentially improv- 
ing the chances of forming planetesimals there (Youdin 
& Chiang, 2004). 

Several recent papers have examined possible modifica- 
tions to radial drift that might arise as a consequence of 
turbulence within the disk (Durisen et al., 2005; Haghigh- 
ipour & Boss, 2003; Rice et al., 2004). The inward mo- 
tion of solid bodies embedded within the disk occurs as 
a consequence of a gas pressure gradient that leads to 
sub-Keplerian gas orbital velocities. In general, radial 
drift drives particles toward pressure maxima, so that 
the motion is inevitably inward in a quiescent disk. In 
a turbulent disk, on the other hand, it may be possible 
to create local pressure maxima that would act as sites 
where solids concentrate. The basic idea is illustrated in 
Figure 19. Such scenarios, although speculative, are po- 
tentially interesting since if the hypothesized local pres- 
sure maximum occurs on a scale Ar, then the local pres- 
sure gradient ~ P/Ar exceeds the global gradient ~ P/r. 
The time scale to concentrate solids locally is then faster 
than the global inspiral time by a factor ~ (Ar/r) 2 . 




FIG. 19 Illustration of how local pressure maxima within a 
disk could concentrate solid bodies, forming a ring in this ide- 
alized axisymmetric example. Local pressure maxima might 
arise as a consequence of turbulence within the disk. 

4. The Goldreich-Ward mechanism 

As we have shown, strong circumstantial evidence sug- 
gests that the formation of planetesimals must occur on 
a time scale that is very short compared to the disk life- 
time. Two hypothesis have been suggested as to how 
that occurs, 

• Planetesimals may form from pairwise collisional 
growth of smaller bodies. In this case, the same 
physical process that allows dust to agglomerate 
into cm sized objects continues uninterrupted up 
to the planetesimal size scale. This is an economi- 
cal hypothesis, and the only difficulty arises due to 
the uncertain sticking efficiency when cm and me- 
ter scale particles collide. Laboratory experiments 
suggest that the probability of particles sticking 
and continuing to grow depends upon details of 
their surface composition (Supulver et al., 1997). 
Guttler et al. (2010) provides a comprehensive re- 
view of the experimental data. 

• Planetesimals may form from the gravitational 
fragmentation of a dense particle sub-disk near 
the midplane of the gas disk. This suggestion - 
made by Goldreich & Ward (1973) 13 - is attractive 
since it forms planetesimals while entirely bypass- 
ing the size scales that are most vulnerable to radial 
drift. However, there are significant theoretical ob- 
jections which render the simplest versions of the 
model untenable. 

In the following, we discuss the Goldreich-Ward (1973) 
mechanism together with recent refinements of the con- 



Similar considerations are discussed in Safronov (1969), who in 
turn quotes earlier work by Gurevich & Lebedinskii from as early 
as 1950. 
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uniformly rotating sheet. 




FIG. 20 Illustration of the Goldreich-Ward mechanism for 
planetesimal formation. A combination of vertical settling 
and (perhaps) radial drift results in a dust sub-disk whose 
density exceeds the local gas density. This sub-disk becomes 
thin enough to be gravitationally unstable, leading to frag- 
mentation into planetesimals. 



cept. The reader should be aware that planetesimal for- 
mation is one of the more uncertain aspects of planet 
formation, and ongoing research continues to study both 
possibilities. 

Figure 20 illustrates the basic idea underlying the 
Goldreich-Ward (1973) mechanism for planetesimal for- 
mation. A combination of vertical settling and radial 
drift of small solid particles results in the formation of 
a dense sub-disk within which the solid density exceeds 
the local gas density (this obviously requires a very thin 
sub-disk if the local ratio of gas to dust surface density 
is comparable to the fiducial global value of 100). The 
solid sub-disk then becomes gravitationally unstable, and 
fragments into bound clumps of solid particles that sub- 
sequently dissipate energy via physical collisions and col- 
lapse to form planetesimals. 



Gravitational instability requires that the disk be mas- 
sive (high surface density) and / or dynamically cold 
(low velocity dispersion). The classic analysis of the 
conditions for gravitational instability is that of Toomrc 
(1964). Here, we consider the stability of a rotating fluid 
sheet - this is somewhat easier than the collisionless cal- 
culation, gives the same answer to a small numerical fac- 
tor when the gas sound speed is identified with the parti- 
cle velocity dispersion, and carries over to the instability 
of a gas disk that we will discuss later. The simplest 
system to analyze is that of a uniformly rotating sheet 
- in what follows I follow the notation and approach of 
Binney & Tremaine (1987). 

The setup for the calculation is as shown in Figure 21. 
We consider a sheet of negligible thickness in the z = 
plane, with constant surface density So and angular 
velocity f2 = ilz. Our aim is to calculate the stability of 
the sheet to in-plane perturbations. Working in a frame 
that corotates with the (unperturbed) angular velocity 
fi, the fluid equations are, 



9S 

at 

~dt 



+ V • (Sv) 
+ (v-V)v 







Vj9 

"IT 



V$ - 2J7 x v 



-n 2 (xx + yy) 



(129) 



(130) 



where the momentum equation picks up terms for the 
Coriolis and centrifugal forces in the rotating frame. 
These equations apply in the z = plane only. The 
gravitational potential $ is given by the Poisson equa- 
tion, 



V 2 $ = 4ttG£(5(z) 



(131) 



which describes $ in all space. In these equations, v = 
v x x + v y y is the velocity in the rotating frame, £ is the 
surface density, and p = p(S) is the vertically integrated 
pressure. The sound speed is denned via, 



dp 
d£ 



(132) 



In the unperturbed state, £ = £ , c E > = < I ) o> v = and 
p = po = p(Y<o). Substituting these values into the mo- 
mentum equation yields V$o = Q 2 (xx + yy). 
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We now consider perturbations to the surface density, 
velocity, pressure and potential, 

E = S + T,i(x,y,t) 

v = V!(x,y,t) 

V = Po+Pi(x,y,t) 

$ = $o + *i(aM/,M) (133) 

where it is assumed that Si <C E etc. Substituting these 
expressions into the fluid equations, and retaining only 
those terms that are linear in the perturbed quantities, 
we find, 



+ EoV-vi =0 



dt 
dvi 

V 2 $x = 4 7 rGS 1 (5(z) 
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-f VSi - V$i - 2fi x vi (135) 
So 
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where we have made use of the fact that since p is only a 
function of E, Vp = (dp/dE)VE. Note that these equa- 
tions only involve temporal or spatial derivatives of the 
perturbed quantities. Since the equations are (by con- 
struction) linear, the evolution of an arbitrary perturba- 
tion can be decomposed into fourier modes. Assuming a 
wavevector k that is parallel to x, we therefore write the 
perturbations in the form, 
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where the final expression describes the potential pertur- 
bations in the z = plane only. Substitution of these 
expressions into the perturbation equations will reduce 
them to algebraic expressions, which can be combined to 
yield the dispersion relation for the system. 

First though, we simplify the system by noting that 
perturbations in E are the source of perturbations in $. 
We can therefore write $ a in terms of E a . To do this, let 
the general form for $i (i.e. not just at z = 0) be, 



i(kx-Lut) 
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where f(z) is some function that needs to be determined. 
Requiring that V 2 <I>i — for z ^ 0, we find, 



d 2 / 
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which has a general solution / = Ae~ kz + B kz , with A 
and B arbitrary constants. Since $i must remain finite 
as z — > ±oo, the general form of $i is, 



= $ e %{kx-wt)-\kz\ _ 

This is valid throughout all space. 
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To determine $ a , we integrate the Poisson equation 
vertically between z = — e and z = +e, 



J + V 2 <5> 1 dz = J + ATrGE 1 S{z)dz. 



(143) 



Mathematically this requires a bit of care, since the inte- 
grand on the left hand side is zero everywhere except at 
z = 0. However, noting that d 2 $i/dx 2 and d 2 $i/dy 2 are 
continuous at z = 0, while d 2 $i/dz 2 is not, we obtain, 
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y 47rGEi<S(z)dz. (144) 
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and, 
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We are now in a position to substitute Ei, vi and $i 
into the remaining equations (continuity plus the x and 
y components of the momentum equation). The resulting 
algebraic equations are, 
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We seek a dispersion relation i.e. a formula for the growth 
rate u = f(k) of modes of different scale k. Eliminating 



v ax and v ay in turn, we obtain, 



c 2 /c 2 -27rGE |A:|±4^ 
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This is the dispersion relation for a uniformly rotating 
thin sheet. The scale-dependence of the different terms 
is shown graphically in Figure 22. 

Looking back to the form of the perturbations, we note 
that the sheet is: 

• STABLE if oj 2 > 0, since in this case u) is real and 
the perturbations are oscillatory. 

• UNSTABLE if cu 2 < 0, for which case uj is imagi- 
nary and perturbations grow exponentially. 

The rotational term (AVt 2 ) is stabilizing at all scales, 
while the pressure term (e 2 fc 2 ) has a strong stabilizing in- 
fluence at large k (i.e. small spatial scales). Self-gravity, 
represented by the — 27rGE |A:| term, has a negative con- 
tribution to oj 2 and so destabilizes the sheet. 

The condition for marginal stability is that cu 2 > at 
all spatial scales. The most unstable scale fc cr it can be 
found by setting dw 2 /dfc = 0, which yields, 
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FIG. 22 The dispersion relation (solid black line) for a uni- 
formly rotating sheet, illustrating the contributions from pres- 
sure, rotation, and self-gravity (dashed blue lines). The sys- 
tem is unstable if, at any value of the wavenumber k, uj 2 falls 
below the red line and is negative. Pressure is a stabiliz- 
ing influence that is most important at large k (small spatial 
scales), while rotation acts to stabilize the system at small k 
(large spatial scales). 



The sheet is marginally stable when w 2 (fc cr i t ) = 0, which 



gives the stability condition as, 
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This analysis can be extended in several ways - for exam- 
ple to include differential rotation or global rather than 
local stability. A generic way of expressing the results of 
such calculations is to define the Toomre Q parameter, 



Q 
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In terms of Q, a disk is unstable 14 to its own self-gravity 
if Q < Q cr it, and stable if Q > Q C rit- Typically Q C rit — 1 
- for the specific system we have investigated it would be 
1/2. 



14 For a differentially rotating disk, it is easy to verify that stability 
depends upon the parameter combination c s Sl/(GSo) via a time 
scale argument. First derive the time scale for shear to separate 
two points that are initially Ar apart, and equate this to the 
collapse time scale under gravity to find the maximum scale on 
which collapse can occur without being affected by shear. Taking 
the ratio of this scale to the Jeans scale (the smallest scale on 
which collapse can occur without being inhibited by pressure 
gradients) yields the correct functional form of Q. 



We have derived the stability of a fluid disk in uniform 
rotation. Differential rotation and global effects alter the 
value of Qcrit, but do not fundamentally change the re- 
sult. For a collisionless disk (e.g. one made of stars or 
small solid particles) a comparable result applies if we re- 
place the sound speed c s by the one-dimensional velocity 
dispersion a. 

The most unstable wavelength is, 



2tt 



2± 



(153) 



Comparing this to the scale height of the disk h — c s /0, 
we find that at marginal stability, 
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i.e. the instability afflicts small-ish spatial scales within 
the disk. 

Let us apply this analysis to the problem of planetes- 
imal formation. If we ignore radial drift, then at 1 AU 
^dust ^ Sgasj or about 10 gem for a minimum 
mass Solar Nebula model (note that a gas to dust ratio 
of 100 is a commonly used approximation in protoplane- 
tary disk theory). Setting Q = ufi/^CJEdust) — L an d 
taking M* = M@ , we find that instability requires a crit- 
ical velocity dispersion in the solid component, 



a ~ 10 cms 



(155) 



Since the gas sound speed at this radius is of the order of 
10 5 cms -1 , and the scale heights of the gas and particle 
disks are respectively proportional to c s and cr, we see 
that an extremely thin disk is required before instability 
will set in! 

If instability occurs, the most unstable wavelength is 
predicted to be, 



Acrit ~ 3 x 10 8 cm. 
The mass within an unstable patch is then, 
to ~ ^SdustA^ - 3 x 10 18 g 
which would correspond to a spherical body of size 
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for a material density of pd = 3 gem 3 . The collapse 
time scale at distance A cr ;t from mass m, 
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is very short - less than a year for the parameters adopted 
above. Even if we allow for the fact that angular momen- 
tum will preclude a prompt collapse, the derived time 
scale for planetesimal formation via gravitational insta- 
bility remains extremely short - perhaps of the order of 
10 3 yr (Goldreich k Ward, 1973). 
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Formation of planetesimals via the Goldreich-Ward 
mechanism has several attractive features, most notably 
the short time scale and complete bypass of the size 
regime most vulnerable to radial drift. However in its 
simplest form, the mechanism fails to work. The problem 
lies in the fact that even in an intrinsically non-turbulent 
gas disk, the formation of a dense solid sub-disk leads 
to self-generated turbulence and associated vertical stir- 
ring prior to gravitational instability. As noted above, 
for gravitational instability to operate we require a thin 
sub-disk in which, for our choice of parameters, 
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Within this midplane layer, the volume density of solids 
would exceed the density of gas by a factor of the order 
of 100 - i.e. the extreme thinness of the solid disk in- 
verts the normal gas to dust ratio which favors gas by 
the same factor. Since the gas and dust are well cou- 
pled for small particle sizes, within the sub-disk (where 
the solid component dominates) we expect both the gas 
and the dust to orbit at the natural velocity for the solid 
component, which is the Kepler velocity. The gas just 
above the layer, on the other hand, will rotate slower due 
to the radial gas pressure gradient. There will therefore 
be a velocity gradient in the z direction that is of the 
order of {h gas /r) Vk //idust- This shear will be Kelvin- 
Helmholtz unstable, leading to turbulence that prevents 
the layer ever getting thin enough to fragment into plan- 
etesimals (Cuzzi, Dobrovolskis & Champney, 1993). The 
condition for Kelvin-Hclmholtz instabilities to develop 
(Sekiya, 1998; Youdin & Shu, 2002) is that the Richard- 
son number, which measures the competition between 
vertical shear and buoyancy, is Ri < Ri cr it, where, 
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and N, the Brunt Vaisala frequency, is defined as, 
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The standard stability analysis obtains a critical Richard- 
son number Ri = 0.25, but both analytic calculations 
including the effect of Coriolis forces, and numerical sim- 
ulations, favor a larger value of around unity (Gomez & 
Ostrikcr, 2005; Johansen, Henning & Klahr, 2006). 



5. Streaming instabilities 

In the classical Goldreich-Ward scenario, turbulence 
within the gas disk (whether intrinsic or self-excited) is 
unequivocally detrimental to the chances of the mech- 
anism working. Recently, however, there has been a 
resurgence of interest in related models which share with 
Goldreich-Ward the central idea that planetesimals form 



via gravitational collapse. The key insight is the real- 
ization that protoplanetary disks may support regions 
within which turbulence acts to locally enhance the ra- 
tio of solids to gas. As a result, patches of the disk may 
become dense enough to collapse into planetesimals (Ga- 
raud & Lin, 2004), even though the global disk has too 
small a dust to gas ratio to support Goldreich-Ward plan- 
etesimal formation. It is useful to distinguish two con- 
ceptually different routes by which disk turbulence might 
assist planetesimal formation: 

1 Passive concentration of particles by tur- 
bulence. It has long been known that particles, 
treated as trace contaminants that are advected 
along with (and drift relative to) the gas, can be- 
come concentrated in turbulent flows. How exactly 
this works depends upon the nature of the turbu- 
lence. Generically, we expect strong fluctuations in 
the particle concentration to occur on scales compa- 
rable to the dissipation scale. Larger scale concen- 
trations are possible if the turbulence is character- 
ized by vortices (Barge & Sommeria, 1995), or if, as 
already noted, the turbulence creates long-lived lo- 
cal pressure maxima within the disk. Cuzzi, Hogan 
& Shariff (2008) has advocated a model which ap- 
peals to strong particle concentration in turbulence 
as a prerequisite for planetesimal formation. 

2 Two-fluid instabilities. These are instabilities 
whose existence relies upon there being two-way 
coupling between the gas and dust fluids. The gas 
does not just passively advect particles, but also 
suffers a back-reaction which becomes significant 
at particle concentrations that are much smaller 
than those needed for the onset of gravitational col- 
lapse. A detailed analysis of the two-fluid system 
shows that it supports linear streaming instabili- 
ties (Bai & Stone, 2010; Goodman & Pindor, 2000; 
Youdin & Goodman, 2005) that can, under circum- 
stances possibly realized in nature, result in strong 
local particle concentration. Streaming instabili- 
ties defy an intuitive physical explanation, but they 
appear to be pervasive provided that the particles 
can attain local densities comparable to those of 
the gas (i.e. provided that the gas disk is quies- 
cent enough). Numerical simulations by Johansen 
et al. (2007) and by Johansen, Youdin & Mac Low 
(2009) suggest that the non-linear evolution of the 
streaming instability can result in gravitational col- 
lapse to form planetesimals, at least provided that 
prior collisional growth manages to form bodies of 
cm-scale. 

A great deal more work needs to be done on all of these 
mechanisms, but the streaming instability, in particular, 
appears to be very promising (for a review see Chiang 
& Youdin, 2009). It may play a central role in plan- 
etesimal formation. Open questions include the extent 
to which large-scale local particle concentration - by ra- 
dial drift (Youdin & Chiang, 2004; Youdin & Shu, 2002) 



.32 



or via gas photoevaporation (Throop & Bally, 2005) - 
may be needed prior to the onset of streaming instabil- 
ities, and how easy it is for collisional evolution to form 
the relatively large particles needed to trigger the insta- 
bility. Coagulation calculations by Zsom et al. (2010), 
for example, have identified a "bouncing barrier" that 
makes growth past mm-scales - which are rather small 
compared to what is needed for streaming instabilities - 
problematic. Empirically, it is worth bearing in mind, 
first, that in the Solar System planetesimals managed 
to form across a wide range of orbital radii, and that, 
second, subtle evidence of the formation channel may 
yet be preserved in the size distribution of asteroids and 
other minor bodies in the Solar System. Morbidelli et al. 
(2009), for example, has argued that the size distribu- 
tion of asteroids favors models in which the initial gener- 
ation of planetesimals were born very large (say 100 km). 
Large planetesimals are the prediction of current models 
based upon the streaming instability, though it is proba- 
bly too early to be entirely confident that this prediction 
is robust. 
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B. Growth beyond planetesimals 

Once planetesimals have formed, further interaction 
between the solid and gaseous components of the disk is 
limited until bodies with sizes > 10 3 km form that are 
large enough to have a gravitational coupling to the gas 15 . 
We will discuss the impact of gravitational coupling ( "mi- 
gration") later in the context of the early evolution of 
planetary systems. If coupling with the gas disk can be 
neglected, further growth to form protoplanets or plan- 
etary embryos is a well-posed N-body problem in which 
gravity provides the dominant physics. 

Being well-posed is not the same as easy - if the Earth 
formed from 5 km radius planetesimals then N ~ 10 9 . 
Although N-body simulations with this many particles 
are certainly feasible (Springel et al., 2005), it is not pos- 
sible to simulate such large particle numbers for the ~ 10 8 
orbits required for an ah initio calculation of terrestrial 
planet formation (making matters more difficult, for long 
duration integrations special numerical techniques are 
needed to keep integration errors under control). The 
usual approach is therefore a combination of statistical 
and N-body methods. 



15 Strictly, all that is known for sure is that aerodynamic effects 
are negligible for planetesimals. The gas disk might still couple to 
planetesimals gravitationally, if turbulence produces surface den- 
sity fluctuations that can gravitationally scatter planetesimals. 
There have been a number of recent studies of this process (John- 
son, Goodman & Menou, 2006; Laughlin, Steinacker & Adams, 
2004; Nelson, 2005; Nelson & Papaloizou, 2004; Yang, Mac Low 
& Menou, 2009). 
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FIG. 23 Setup for calculation of gravitational focusing. Two 
bodies of mass m, moving on a trajectory with impact pa- 
rameter b, have a velocity at infinity of a/2. 



1. Gravitational focusing 

For sufficiently small bodies, the effects of gravity can 
be ignored for the purposes of determining whether they 
will physically collide. A massive planet, on the other 
hand, can gravitationally focus other bodies toward it, 
and as a result has a collision cross section that is much 
larger than its physical cross section. 

To evaluate the magnitude of this gravitational focus- 
ing, consider two bodies of mass m, moving on a trajec- 
tory with impact parameter b, as shown in Figure 23. The 
relative velocity at infinity is a. At closest approach, the 
bodies have separation R c and velocity V max . Equating 
energy in the initial (widely separated) and final (closest 
approach) states we have, 
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Noting that there is no radial component to the veloc- 
ity at the point of closest approach, angular momentum 
conservation gives, 
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If the sum of the physical radii of the bodies is R s , then 
for R c < R s there will be a physical collision, while larger 
R c will result in a harmless flyby 16 . The largest value of 
the impact parameter that will lead to a physical collision 
is thus, 



b z = Ri 



AGmR s 
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which can be expressed in terms of the escape velocity 
from the point of contact, v 2 sc = 4Gm/R s as, 



R 1 
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16 This is true for solid bodies — for giant planets or stars tidal 
effects can lead to significant dissipation of energy even when 
R c > R s (Fabian, Pringle & Rees, 1975). 
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The cross section for collisions is then, 

T = irR 2 s (l + ^) , (167) 

where the term in brackets represents the enhancement to 
the physical cross section due to gravitational focusing. 
Clearly a planet growing in a "cold" planetesimal disk 
for which a <C v csc will grow much more rapidly as a 
consequence of gravitational focusing. As a consequence, 
determining the velocity dispersion of bodies of different 
masses during the planet formation process is extremely 
important. 

2. Growth versus fragmentation 

When two initially solid bodies physically collide the 
outcome can be divided broadly into three categories: 

• Accretion. All or most of the mass of the impactor 
becomes part of the mass of the final body, which 
remains solid. Small fragments may be ejected, but 
overall there is net growth. 

• Shattering. The impact breaks up the target 
body into a number of pieces, but these pieces re- 
main part of a single body (perhaps after reaccumu- 
lating gravitationally). The structure of the shat- 
tered object resembles that of a rubble pile. 

• Dispersal. The impact fragments the target into 
two or more pieces that do not remain bound. 

To delineate the boundaries between these regimes quan- 
titatively, we consider an impactor of mass m colliding 
with a larger body of mass M at velocity v. We define 
the specific energy Q of the impact via, 



and postulate, plausibly, that this parameter largely con- 
trols the result. The thresholds for the various collision 
outcomes can then be expressed in terms of Q. Conven- 
tionally, we define the threshold for catastrophic disrup- 
tion Q* D as the minimum specific energy needed to dis- 
perse the target in two or more pieces, with the largest 
one having a mass M/2. Similarly Q* s is the threshold 
for shattering the body. More work is required to dis- 
perse a body than to shatter it, so evidently Q* D > Q* 5 . 
It is worth keeping in mind that in detail the outcome 
of a particular collision will depend upon many factors, 
including the mass ratio between the target and the im- 
pactor, the angle of impact, and the shape and rotation 
rate of the bodies involved. Quoted values of Q* D are of- 
ten averaged over impact angles, but even when this is 
done the parameterization of collision outcomes in terms 
of Q is only an approximation. 

The estimated values of Q * D for a target of a particular 
size vary by more than an order of magnitude depending 



upon the composition of the body, which can broadly 
be categorized into solid or shattered rock, and solid or 
porous ice. For any particular type of body, however, two 
distinct regimes can be identified: 

• Strength dominated regime. The ability of 
small bodies to withstand impact without being 
disrupted depends upon the material strength of 
the object. In general, the material strength of 
bodies declines with increasing size, owing to the 
greater prevalence of defects that lead to cracks. In 
the strength dominated regime Q* D decreases with 
increasing size. 

• Gravity dominated regime. Large bodies are 
held together primarily by gravitational forces. In 
this regime Q* D must at the very least exceed the 
specific binding energy of the target, which scales 
with mass M and radius a as Qb oc GM/a oc p^a 2 . 
In practice it requires a great deal more than this 
minimum amount of energy to disrupt the target - 
so Qb is not a good estimate of Q * D - but nonethe- 
less Q* D does increase with increasing size. 

Although the transition between these regimes is reason- 
ably sharp there is some influence of the material proper- 
ties (in particular the shear strength) on the catastrophic 
disruption threshold for smaller bodies within the gravity 
dominated regime. 

Values of Q* s and Q* D can be determined experimen- 
tally for small targets (Arakawa, Leliwa-Kopystynski & 
Maeno, 2002). Experiments are not possible in the grav- 
ity dominated regime, but Q* D can be estimated theo- 
retically using numerical hydrodynamics (Benz & As- 
phaug, 1999; Leinhardt & Stewart, 2009) or (for rub- 
ble piles) rigid body dynamics simulations (Korycansky 
k Asphaug, 2006; Leinhardt & Richardson, 2002). The 
simplest parameterization of the numerical results is as 
a broken power law that includes terms representing the 
strength and gravity regimes, 

% = ^(l^) C + ^(l^) d - (169) 

Often (but not always) Q* D is averaged over impact ge- 
ometry, and q s , q g , c and d are all constants whose values 
are derived by fitting to the results of numerical simula- 
tions. 

Benz & Asphaug (1999) and Leinhardt & Stewart 
(2009) determined the values of the fitting parameters in 
equation (169) from the results of an ensemble of simu- 
lations of impacts into icy or rocky targets. Their results 
are given in Table III and plotted as a function of target 
size in Figure 24. One observes immediately that the re- 
sults for a particular target material vary with the impact 
velocity, and hence that Q* D is not the sole determinant 
of the outcome of collisions. There is, however, a clear 
transition between the strength and gravity dominated 
regimes, with the weakest bodies being those whose size 
is comparable to the cross-over point. The most vulner- 
able bodies are generally those with radii in the 100 m to 
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TABLE III Parameters for the catastrophic disruption 
threshold fitting formula (equation 169), which describes how 
Q* D scales with the size of the target body. The quoted val- 
ues were derived by Benz & Asphaug (1999) and Leinhardt & 
Stewart (2009) using numerical hydrodynamics simulations of 
collisions, which are supplemented in the strength dominated 
regime by experimental results. 
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1 km range. Just how vulnerable such bodies are to catas- 
trophic disruption depends sensitively on their make-up, 
and it would be unwise to place too much trust in pre- 
cise numbers. As a rough guide, however, the weakest 
icy bodies have minimum Q* D ~ 10 5 erg g _1 , while the 
strongest conceivable planetesimals (unfractured rocky 
bodies) have minimum Q* D > 10 6 erg g _1 . 

As a reality check, we may note that asteroids in the 
main belt with e ~ 0.1 would be expected to collide to- 
day with typical velocities of the order of 2 km s _1 . For a 
mass ratio m/M = 0.1 the specific energy of the collision 
is then around Q = 2 x 10 9 erg g _1 , which from Figure 24 
is sufficient to destroy even quite large solid bodies with 
a ~ 100 km. This is consistent with the observation of as- 
teroid families, and the interpretation of such families as 
collisional debris. Evidently the random velocities that 
characterize collisions must have been much smaller dur- 
ing the epoch of planet formation if we are to successfully 
build large planets out of initially km-scale planetesimals. 



3. Shear versus dispersion dominated encounters 

A more subtle distinction that nevertheless plays a 
crucial role in planet formation is whether encounters 
between bodies can be described via 2-body dynamics 
- in which only the gravity of the two objects them- 
selves matters — or whether the tidal influence of the 
Sun also needs to be considered (3-body dynamics) . Gol- 
dreich, Lithwick & Sari (2004) have recently summarized 
in simple terms why the distinction between 2 and 3- 
body dynamics matters at different stages of the planet 
formation process. We consider a 3-body system consist- 
ing of a large body (a "planet") with mass M, a small 
body of negligible mass (described as a test particle), 
and the Sun, and define the Hill radius r# as the radius 
within which the gravity of the planet dominates (in as- 
trophysical contexts, the same concept is referred to as 
the "Roche lobe"). Roughly, this is obtained by equating 
the angular velocity for an orbit at distance rn from the 
planet with the angular velocity of the planet around the 
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FIG. 24 The specific energy Q* D for catastrophic disruption of 
solid bodies is plotted as a function of the body's radius. The 
solid and short dashed curves show results obtained using fits 
to theoretical calculations for impacts into "strong" targets 
by Benz & Asphaug (1999). The long dashed curve shows 
the recommended curve for impacts into "weak" targets from 
Leinhardt & Stewart (2008), derived from a combination of 
impact experiments and numerical simulations. In detail the 
solid curves show results for basalt at impact velocities of 
5 km s _1 (upper curve) and 3 km s -1 (lower curve). The short 
dashed curves show results for water ice at 3 km s _1 (the lower 
curve for small target sizes) and 0.5 km s _1 (upper curve for 
small target sizes). The long dashed curve shows results for 
normal impacts into weak water ice targets at 1 km s _1 . 

star. We find, 

where the factor 3 is included for consistency with more 
detailed derivations. For circular orbits, collisions are 
forbidden for an orbital separation Ao between the small 
body and the planet such that Aa < r# (c.f. the Trojan 
asteroids in the Solar System). If we define a character- 
istic velocity at the Hill radius, 



then for, 

• a > v H 2-body dynamics describes collisions quite 
well. This regime is called dispersion domi- 
nated. 

• u < Vh 3-body effects are important. This regime 
is called shear dominated. 
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When a < vh and we are shear dominated, the collision 
rate is modified compared to expectations based on 2- 
body dynamics. 



The radius of the planet grows at a linear rate. Writing 
the planet mass M = (4/3)7i\Rj}p planet) where p p i anc t is 
the planet density, 



4. Growth rates 
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We now proceed to derive an estimate for how fast a 
planet will grow due to accretion of planetesimals. We 
assume that the growing body, of mass M, radius R s , and 
surface escape speed v esc is embedded within a "swarm" 
of planetesimals with local surface density £ p , velocity 
dispersion a, and scale height h p . The volume density of 
the planetesimal swarm is, 



2K 



(172) 



Then if 3-body effects can be ignored, the large body 
grows at a rate, 



—^j- = PswCnrH s 11 + -^2- 



(173) 



This can be simplified since hp ~ o~/£l and hence p s 
inversely proportional to a. We find, 



is 



1 + s?) (174) 



dM 1 n 2 



where the numerical prefactor, which has not been de- 
rived accurately here, depends upon the assumed veloc- 
ity distribution of the planetesimals. For an isotropic 
distribution the prefactor is \/3/2 (Lissauer, 1993). 
This simple result is important. We note that: 

• The velocities of the planetesimals enter only via 
the gravitational focusing term, which can however 
by very large. 

• The rate of mass growth scales linearly with £ p — 
we expect faster growth in disks that have more 
mass in planetesimals (due to a higher gas mass 
and / or a higher ratio of solids to gas). 

• Other things being equal, growth will be slower at 
large radii, due to lower E p and smaller fi. 

Complexity arises because as a planet grows, it stars to 
influence both the velocity dispersion and, eventually, the 
surface density of the planetesimal swarm in its vicinity. 

Two simple solutions of the growth equation give an 
idea of the possibilities present in more sophisticated 
models. First, assume that the gravitational focusing 
term 



F g is constant. In this regime, 

^ oc i? 2 oc MV3 
dt s 



which has solution, 



R s oc t. 



(175) 



(176) 



If we assume that at the orbital radius of Jupiter £ p 
10 gem -2 , then for p p i anct = 3 gem -3 , 



~ 0.2F„ cm yr . 

dt 9 



(178) 



This initial growth rate is slow, which implies that to 
form the cores of the giant planets in a reasonable time, 
large gravitational focusing factors are needed. For ex- 
ample, to reach 1000 km in 10 5 yr, we require F g <~ 5000. 
The need for large gravitational enhancements to the col- 
lision rate is even more severe for the ice giants, but sub- 
stantially easier in the terrestrial planet region. 

Since empirically F g must be large, a second useful 
limit to consider is the case where F g ^> 1. If we as- 
sume that a is constant (i.e. consider the regime where 
the growing planet has not yet managed to dominate the 
dynamical excitation of the planetesimal swarm) then, 
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The growth equation (174) gives, 



dM 
dt 



oc MR, 



with solution, 



M 



(M " 1/3 - ktf' 



(179) 



(180) 



(181) 



where Mq is the initial mass at time t = and k is 
a constant. In this regime the increasing gravitational 
focusing factor means that M — > oo in a finite time, 
allowing much more rapid growth. 



5. Isolation mass 

As noted above, rapid growth requires that a remain 
low — i.e. that the planetesimals remain on roughly 
circular orbits. This means that there is a finite supply 
of planetesimals that have orbits that pass close enough 
to a growing planet to collide — once these have all been 
consumed growth is bound to slow. The mass at which 
this slowdown occurs is described as the isolation mass 
M iso . 
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To estimate the isolation mass, we note that a planet 
grows by accreting planetcsimals within a 'feeding zone'. 
The size of the feeding zone Aa max is set by the maximum 
distance over which the planet's gravity is able to perturb 
planetesimal orbits sufficiently to allow collisions, so it 
will scale with the Hill radius. Writing 



Cr 



H 



(182) 



with C a constant of order unity, we have that the mass 
of planetesimals within the feeding zone is, 



2na ■ 2Aa n 



(183) 



Note the 1/3 power of the planet mass, which arises 
from the mass dependence of the Hill radius. As a 
planet grows, its feeding zone expands, but the mass of 
new planetesimals within the expanded feeding zone rises 
more slowly than linearly. We thus obtain the isolation 
mass by setting the planet mass equal to the mass of the 
planetesimals in the feeding zone of the original disk, 



M iso = 47ra • C 
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which gives, 



Miso = ±^/*C 3 ' 2 M^ 1/2 ?*/ 2 a 3 . 
v3 



(184) 



(185) 



Evaluating this expression in the terrestrial planet region, 
taking a = 1 AU, S p = 10 gem" 2 , M* = M Q and C = 
2^3 (L issauer, 1993), we obtain, 



M iso ~ 0.07 M a 



(186) 



Isolation is therefore likely to occur late in the forma- 
tion of the terrestrial planets. Repeating the estimate for 
the conditions appropriate to the formation of Jupiter's 
S p = 10 gcm~ 2 as adopted by Pollack et al. 
gives, 



(1996) 
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M iso ~ 9 M a 



(187) 



This estimate is comparable to, or larger than, the cur- 
rent best determinations for the mass of the Jovian core 
(Guillot, 2005). Full isolation may or may not be relevant 
to the formation of Jupiter, depending upon the adopted 
disk model. 



6. Coagulation equation 

One might legitimately question whether the assump- 
tion that the mass distribution of growing bodies can 



17 Note that this is a factor of several enhanced above the minimum 
mass Solar Nebula value. 



be neatly divided into two groups — planetesimals and 
growing planetary embryos — is any good. The quan- 
titative approach to describing the evolution of a arbi- 
trary size distribution is based on the coagulation equa- 
tion (Smoluchowski, 1916). This allows us to drop the 
two groups approximation though at the expense of an 
enormous increase in complexity. 

To write the coagulation equation in its simplest 
form 18 , assume that the masses of bodies are integral 
multiples of some small mass mi. At time t there are n k 
bodies of mass m,k — km\. The coagulation equation in 
discrete form is, 



dt ~ 2 Is A ^ 



J2 A k 



(188) 



i+j=k 



where Aij is the rate of mergers between bodies of mass 
rrii and nij. The first term on the right-hand side of the 
equation describes the increase in the number of bodies 
of mass nrik due to collisions of all possible pairs of bodies 
whose masses m^ and rrij sum to m,k- The second term 
describes the decrease due to bodies of mass nik being 
incorporated into even larger bodies. The possibility of 
fragmentation is here neglected. In this formulation of 
the problem of planetary growth, all of the physics - 
such as gravitational focusing — enters via the rate co- 
efficients A^. 

Equation (188), or variants of it, has been used ex- 
tensively to study planet formation (Inaba et al., 2001; 
Kcnyon & Luu, 1998; Safronov, 1969; Wetherill & Stew- 
art, 1993), either on its own or in combination with direct 
N-body simulations (Bromley & Kcnyon, 2006). Gener- 
ally the coagulation equation needs to be supplemented 
with additional equations that describe the evolution of 
the velocity dispersion as a function of mass, as described 
for example in Kenyon & Luu (1998). Because of the fact 
that all i, j such that m t + m,j = mk contribute to the 
evolution of n& , even the coagulation equation on its own 
is not a simple equation to deal with, and few analytic 
solutions are known. One (over)-simple case for which 
an analytic solution exists is for the case when, 

A^ = a (189) 

with a a constant. Then, if the initial state of the sys- 
tem comprises n\ bodies all of mass mi, the solution to 
equation (188) is, 

n k = m/ 2 (l-/) fe - 1 
/ = 1 • (190) 



18 It is also possible to write the coagulation equation as an integro- 
differcntial equation for a continuous mass function n(m, t) 
(Safronov, 1969), or as a discrete equation where bodies are 
binned into arbitrary mass intervals (typically logarithmic). 
Kenyon & Luu (1998) provide a clear description of how the 
coagulation equation may be formulated and solved in the more 
general case. 
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FIG. 25 Illustrative analytic solution to the coagulation equa- 
tion for the simple case in which Aij = a, with a a constant. 
Initially all bodies have mass mi. The solution is plotted for 
scaled times t' = anit equal to 1, 10, 100 and 10 . The up- 
per panel shows the number of bodies n k of each mass (the 
vertical scale is arbitrary), while the lower panel shows how 
the mass distribution evolves. This solution is an example of 
orderly growth — as time progresses the mean mass steadily 
increases while the shape of the mass spectrum remains fixed. 

This solution is shown as Figure 25. The mass spectrum 
remains smooth and well-behaved as growth proceeds, 
and with increasing time the characteristic mass increases 
linearly while maintaining a fixed shape. 

More generally, solutions to the coagulation equation 
fall into two classes (e.g. Lee, 2000): 

• Solutions that exhibit orderly growth, in which the 
mass distribution evolves smoothly with time to- 
ward higher mean masses. The analytic solution 
given above for the case = constant is an ex- 
ample of this type of evolution. Another analytic 
example is Aij cx (rrij + mj ) . 

• Solutions that show runaway growth. In this case 
the mass distribution develops a power-law tail to- 
ward high masses — physically this corresponds to 
one or a handful of bodies growing rapidly at the 
expense of all the others. The long-term validity 
of the coagulation equation once runaway growth 
occurs is evidently limited. An analytic example 
occurs for a rate coefficient A^ cx mirrij. 

Looking back to equation (174), we note that the rate 
coefficient is expected to scale as A cx i? 2 cx to 2 / 3 in the 
regime where gravitational focusing is unimportant, and 
A cx i? 2 v 2 sc cx m 4 / 3 once gravitational focusing is dom- 
inant. By comparison with the aforementioned analytic 



solutions, we expect that the initial growth of planetes- 
imals will occur in the orderly regime, while runaway 
growth may occur once the largest bodies are massive 
enough for gravitational focusing to become significant. 



7. Overview of terrestrial planet formation 

We conclude the discussion of terrestrial planet for- 
mation by summarizing briefly the main stages of the 
process: 

1. Dust particles agglomerate to form, eventually, 
planetesimals. Initially this almost certainly occurs 
via pairwise collisions, though how (or whether) 
this process can continues to work for cm to meter 
scale particles remains somewhat murky. Gravita- 
tional instability may allow a bypass of these tricky 
sizes. 

2. Growth beyond planetesimals occurs via direct col- 
lisions, with an increasing role for gravitational fo- 
cusing as masses become larger. Dynamical friction 
keeps the velocity dispersion of the most massive 
bodies low. A phase of runaway growth occurs in 
which a few bodies grow rapidly at the expense of 
the rest. 

3. Runaway growth ceases once the largest bodies be- 
come massive enough to stir up the planetesimals in 
their vicinity. A phase of oligarchic growth ensues, 
in which the largest objects grow more slowly than 
they did during runaway growth, but still more 
rapidly than small bodies (Kokubo & Ida, 1998; 
Thommes, Duncan & Levison, 2003). Growth con- 
tinues in this mode until the isolation mass is ap- 
proached, at which point growth slows further. 

4. Further evolution occurs as a result of collisions 
between the initially relatively isolated planetary 
embryos left over after oligarchic growth. The em- 
bryos are perturbed onto crossing orbits due to the 
influence of the giant planets and mutual secular 
resonances (Chambers & Wetherill, 1998). The fi- 
nal assembly of the terrestrial planets takes around 
100 Myr, with the predicted configuration vary- 
ing depending upon the assumed surface density of 
planetesimals and existence (or not) of giant plan- 
ets (Kokubo, Kominami & Ida, 2006; Levison & 
Agnor, 2003; Raymond, Quinn & Lunine, 2005). In 
the Solar System, one of the final impacts on the 
Earth is widely considered to have given rise to the 
ejection of enough mass into orbit to subsequently 
form the Moon (Canup, 2004). 

The dominant uncertainties in theoretical models for ter- 
restrial planet formation are arguably found during stage 
1 — the formation of planetesimals. It is also true that 
most simulations of the late stages of terrestrial planet 
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formation lead to planetary properties (such as the ec- 
centricity, and the mass of Mars compared to the other 
terrestrial planets) that differ somewhat from those ob- 
served in the Solar System. Thus, although there is gen- 
eral confidence that the basic physics of terrestrial planet 
formation is understood, it is clear that current models 
do not include all of the ingredients needed to accurately 
match Solar System constraints (Raymond et al., 2009). 



C. Gas giant formation 

Two theoretical models vie to explain the formation 
of gas giant planets. The core accretion model (Boden- 
heimer & Pollack, f986; Mizuno, f980), developed in its 
most refined form by Pollack et al. (1996), postulates that 
the envelopes of gas giants are accreted subsequent to the 
formation of a large core, which is itself assembled in a 
manner analogous to terrestrial planet formation. Core 
accretion is the dominant theory for massive planet for- 
mation. The gravitational instability model, on the other 
hand, is based on the idea that a massive protoplane- 
tary disk might collapse directly to form massive planets 
(Cameron, 1978; Kuiper, 1951). Boss (1997) is the most 
recent advocate of this idea, which has come under re- 
newed theoretical scrutiny with the discovery of many 
extrasolar planets with masses much larger than that of 
Jupiter. 

In this Section, we review the physics of these theories 
in turn. We also discuss the observational constraints on 
the different theories, which include inferences as to the 
core masses of the gas giants in the Solar System, the host 
metallicity / planet frequency correlation for extrasolar 
planetary systems, and — indirectly — comparison of 
the theoretically derived time scales with observations 
of protoplanetary disk lifetimes. This is a critical issue, 
since gas giants must form prior to the dispersal of the gas 
disk. Any successful model of massive planet formation 
must grow such bodies within at most 5-10 Myr (Haisch, 
Lada & Lada, 2001). 



1. Core accretion model 

The main stages in the formation of a gas giant via 
core accretion are illustrated in Figure 26. A core of 
rock and / or ice forms via the same mechanisms that we 
have previously outlined for terrestrial planet formation. 
Initially, there is either no atmosphere at all (because 
the potential is too shallow to hold on to a bound at- 
mosphere), or any gas is dynamically insignificant. How- 
ever, as the core grows, eventually it becomes massive 
enough to hold on to a significant envelope. At first, 
the envelope is able to maintain hydrostatic equilibrium. 
The core continues to grow via accretion of planetesi- 
mals, and the gravitational potential energy liberated as 
these planetesimals rain down on the core provides the 
main source of luminosity. This growth continues until 




FIG. 26 Illustration of the main stages of the core accretion 
model for giant planet formation. 

the core reaches a critical mass. Once the critical mass 
is reached, the envelope can no longer be maintained in 
hydrostatic equilibrium. The envelope contracts on its 
own Kelvin-Helmholtz time scale, and a phase of rapid 
gas accretion occurs. This process continues until (a) the 
planet becomes massive enough to open up a gap in the 
protoplanetary disk, thereby slowing down the rate of gas 
supply, or (b) the gas disk itself is dispersed. 

The novel aspect of the core accretion model is the 
existence of a critical core mass. Mizuno (1980) used nu- 
merical models to demonstrate the existence of a maxi- 
mum core mass, and showed that it depends only weakly 
on the local properties of the gas within the protoplane- 
tary disk. A clear exposition of this type of calculation 
is given in, for example, Papaloizou & Terquem (1999). 
Here, following Stevenson (1982), we show that a toy 
model in which energy transport is due solely to radia- 
tive diffusion displays the key property of a critical core 
mass. 

Consider a core of mass M cole and radius i? C ore, sur- 
rounded by a gaseous envelope of mass M env . The total 
mass of the planet, 

M t = M core + M env . (191) 

The envelope extends from i? coro to some outer radius 
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i? out , which marks the boundary between the gas bound 
to the planet and the gas in the protoplanetary disk. 
i?out may be determined by thermal effects (in which 
case i?out ~ GM t /c 2 s , with c s the disk sound speed) or 
by tidal considerations (giving an outer radius of m), 
whichever is the smaller. If the envelope is of low mass, 
then the largest contribution to the luminosity is from 
accretion of planetesimals onto the core. This yields a 
luminosity, 



L = 



GM nore M, 



corc lvJ corc 



R c 



(192) 



which is constant through the envelope. 

If we assume that radiative diffusion dominates the 
energy transport, then the structure of the envelope is 
determined by the equations of hydrostatic equilibrium 
and radiative diffusion, 
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Having derived the density profile the mass of the enve- 
lope follows immediately, 
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The right-hand-side of this equation has a strong depen- 
dence on the total planet mass M t and a weaker depen- 
dence on the core mass M coro via the expression for the 
luminosity, 



L = GMc „° rcMcorc oc MHlMc, 



1 core 1 "core- 



(202) 



In principle there are further dependencies to consider 
since -R ou t is a function of M t and R CO re is a function of 
Mcorc, but these enter only logarithmically and can be 
safely ignored. Noting that, 



M rnrr = M t — M e 



(203) 



we find that, 



where a is the Stcfan-Boltzmann constant and kr the 
Rosseland mean opacity (assumed constant). Eliminat- 
ing the density between these equations we find that, 
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We now integrate this equation inward from the outer 
boundary, making the approximation that M(r) w M t 
and taking L and kr to be constants, 



r T 3 dT = 3kr J; r dP . 



(196) 



Once we are well inside the planet we expect that T A 3> 
T^ isk and that P 3> Pdisk, so the integral yields, approx- 
imately, 
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Substituting P in this equation with an ideal gas equation 
of state, 



firrip 



(198) 



we eliminate T 3 in favor of the expression involving 
dT/dr and integrate once more with respect to radius 
to obtain, 
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where we have shown explicitly the dependence on the en- 
velope opacity and planetesimal accretion rate but have 
swept all the remaining constants (and near-constants) 
into a single constant C. 

Solutions to equation (204) are plotted as Figure 27. 
One sees that for fixed M core , there exists a maximum or 
critical core mass M cr i t beyond which no solution is pos- 
sible. The physical interpretation of this result — whose 
origin is not terribly clear even within this simple model 
— is that if one tries to build a planet with a core mass 
above the critical mass hydrostatic equilibrium cannot be 
achieved in the envelope. Rather the envelope will con- 
tract, and further gas will fall in as fast as gravitational 
potential energy can be radiated. 

This toy model should not be taken too seriously. How- 
ever, it does illustrate the most important result from 
more detailed calculations — namely that the critical 
mass increases with larger M coro and with enhanced opac- 
ity. Ikoma, Nakazawa & Emori (2000) derive an approx- 
imate fit to numerical results, 
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where the power-law indices are uncertain by around 
±0.05. The weak dependence of the critical core mass on 
the planetesimal accretion rate means that, within a par- 
ticular core accretion model, we can always speed up the 
approach to runaway gas accretion simply by increasing 
the assumed surface density of planetesimals in the vicin- 
ity of the growing core. Contrary to what is sometimes 
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FIG. 27 Solutions to equation (204) for the core mass M core 
and total mass M to tai- The blue curve is for a higher plan- 
etesimal accretion rate than for the red curve. The critical 
core mass is shown as the vertical dashed line. One should 
not take solutions to this toy model very seriously, but the 
numbers have been fixed here to correspond roughly to the 
values obtained from real calculations. 



implied, there is no intrinsic difficulty in building planets 
quickly via core accretion. However, faster growth oc- 
curs at the expense of a larger final core mass. As we will 
shortly note, this tradeoff is of concern since estimates 
of the core mass of Jupiter are smaller than the values 
obtained in the classic calculations of core accretion by 
Pollack et al. (1996). 

Although they appear very detailed, extant calcula- 
tions of planet growth via core accretion should probably 
be regarded as illustrative rather than definitive. Two 
sources of uncertainty are particularly worrying: 

• What is the magnitude of the opacity? Al- 
though kr enters equation (205) as rather a weak 
power, its magnitude is highly uncertain. Hubickyj, 
Bodcnhcimer & Lissauer (2005), and more recently 
Movshovitz et al. (2010), have computed core ac- 
cretion models in which the opacity is either arbi- 
trarily reduced or computed from a settling and co- 
agulation model. These models suggest, first, that 
the appropriate value of the opacity in the enve- 
lope is greatly reduced (by a factor of the order 
of 10 2 ) from the interstellar value (Podolak, 2003). 
Second, they indicate that the reduced opacity re- 
sults in substantially faster growth of massive plan- 
ets. Formation time scales as short as a Myr, or 
(for longer formation times) core masses as small 
as 5 Mg, now appear achievable. 



• The neglect of core migration. Theoretical 
work, which we will discuss more fully in a sub- 
sequent Section, suggests that planets or planetary 
cores with masses exceeding M§ are highly vulner- 
able to radial migration as a consequence of grav- 
itational torques exerted by the gas disk. This ef- 
fect is not included in the calculations of Pollack 
et al. (1996) or Hubickyj, Bodenheimer & Lissauer 
(2005). Papaloizou & Terquem (1999) and Alibert 
et al. (2005) have studied the effect of steady in- 
ward migration on core formation, and have showed 
that it makes a large change to the time scale and 
outcome of the process. Matters could be different 
again if the migration process is instead unsteady 
(Rice & Armitage, 2003), or outward. Radial mi- 
gration could also be driven by dynamical inter- 
actions between growing cores and planetesimals 
(Levison, Thommes & Duncan, 2010). 

To summarize, the broad outlines of how core accretion 
works are well established, but even the most sophisti- 
cated models are probably lacking some essential physical 
elements. 



2. Gravitational instability model 

A sufficiently massive and / or cold gas disk is gravita- 
tionally unstable 19 . If — and this is a big if — gravita- 
tional instability leads to fragmentation this can lead to 
massive planet formation (Cameron, 1978; Kuiper, 1951). 
Durisen et al. (2007) provides a recent review of the sta- 
tus of the gravitational instability model for giant planet 
formation. 

We have already derived the necessary conditions for 
gravitational instability to occur. We need the Toomre 
Q parameter to be low enough, specifically, 

Q = ^ < Qcrit - 1 (206) 

where c s is the sound speed in a gas disk of local sur- 
face density S and the disk mass is assumed to be 
small enough that the distinction between the orbital 
and epicyclic frequencies is of little import. If we con- 
sider a disk with h/r — 0.05 at 10 AU around a Solar 
mass star, then the relation h/r = c s /v$ yields a sound 
speed c s ~ 0.5 kms -1 . To attain Q = 1, we then require 
a surface density, 

S w 1.5 x 10 3 gem 2 . (207) 



The terminology used to discuss this process is potentially con- 
fusing. I will use the term gravitational instability to refer to 
disks in which the self-gravity of the gas is significant enough to 
alter the structure or evolution of the disk. Fragmentation refers 
to the case where gravitational instability leads to the breakup 
of the disk into bound objects. 
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This is much larger than estimates based, for example, 
on the minimum mass Solar Nebula, from which we con- 
clude robustly that gravitational instability is most likely 
to occur at an early epoch when the disk mass is still high. 
Recalling that the characteristic wavelength for gravita- 
tional instability is A cr it = 2c 2 / (G£), we find that the 
mass of objects formed if such a disk fragmented would 
be, 
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where Mj is the mass of Jupiter. These order of magni- 
tude estimates suffice to indicate that gravitational insta- 
bility followed by fragmentation could form gas giants. 

It is also straightforward to derive where in the disk 
gravitational instability is most likely to occur. Noting 
that in a steady-state accretion disk z/£ = M/(3n), we 
use the a prescription (Shakura & Sunyaev, 1973) and 
obtain, 



M 



(209) 



The sound speed in a protoplanetary disk decreases 
outward, so a steady-state disk becomes less stable at 
large radii. Indeed, unless the temperature becomes so 
low that external irradiation (not that from the central 
star) dominates the heating, a steady-state disk will be- 
come gravitational unstable provided only that it is large 
enough. 

To derive sufficient conditions for fragmentation, we 
need to go beyond these elementary considerations and 
ask what happens to a massive disk as instability is ap- 
proached. The critical point is that as Q is reduced, non- 
axisymmetric instabilities set in which do not necessarily 
lead to fragmentation. Rather, the instabilities gener- 
ate spiral arms (Laughlin & Bodenheimer, 1994) which 
both transport angular momentum and lead to dissipa- 
tion and heating. The dissipation in particular results 
in heating of the disk, which raises the sound speed and 
makes fragmentation less likely. On a longer time scale, 
angular momentum transport also leads to lower surface 
density and, again, enhanced stability (Lin & Pringle, 
1990). 

Given these consideration, when will a disk fragment? 
Gammie (2001) used both analytic arguments and local 
numerical simulations to identify the cooling time as the 
control parameter determining whether a gravitationally 
unstable disk will fragment. For an annulus of the disk 
we can define the equivalent of the Kelvin- Helmholtz time 
scale for a star, 
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where U is the thermal energy content of the disk per 
unit surface area. Then for an ideal gas equation of state 
with 7 = 5/3 the boundary for fragmentation is: 

• icooi ^ 3S1 1 — the disk fragments. 



• t coo \ > 30 — 1 — disk reaches a steady state in which 
heating due to dissipation of gravitational turbu- 
lence balances cooling. 

This condition is intuitively reasonable. Spiral arms re- 
sulting from disk self-gravity compress patches of gas 
within the disk on a time scale that is to order of mag- 
nitude il^ 1 . If cooling occurs on a time scale that is 
shorter that O -1 , the heating due to adiabatic compres- 
sion can be radiated away, and in the absence of extra 
pressure collapse is likely. It is also worth noting that 
although the above condition was derived locally, global 
simulations show that it provides a good approximation 
to the stability of protoplanetary disks more generally 
(Rice et al., 2003b). One can also express the fragmenta- 
tion boundary in terms of a maximum stress that a self- 
gravitating disk can sustain without fragmenting (Gam- 
mie, 2001). Writing this in terms of an effective a pa- 
rameter, a max ~ 0.06 (Rice, Lodato & Armitage, 2005). 

In a real disk, the cooling time is determined by the 
opacity and by the mechanism of vertical energy trans- 
port: either radiative diffusion or convection. Using 
a disk model, one can then estimate analytically the 
conditions under which a disk will become unstable to 
fragmentation (Clarke, 2009; Levin, 2007; Rafikov, 2005, 
2009). For standard opacities, the result is that fragmen- 
tation is expected only at quite large radii of the order 
of 50 or 100 AU. At smaller radii the disk may still be 
gravitationally unstable, but the instability will saturate 
prior to fragmentation and, instead, contribute to angu- 
lar momentum transport. The upshot is that whether 
fragmentation will occur boils down to a question about 
the size of the disk, which is itself determined by the 
star formation process. Molecular cloud cores that have 
a significant degree of rotational support (J3 > 10~ 2 ) will 
collapse to form disks that are potentially large enough 
to fragment at their outer edges, whereas less rapidly 
rotating cores (or cores that experience significant mag- 
netic braking during collapse) yield compact disks that 
are everywhere stable (Rice, Mayo & Armitage, 2010). 

Most, but not all, numerical simulations yield results 
for disk fragmentation that are broadly consistent with 
the aforementioned analytic arguments 20 . Calculations 
by Boley et al. (2010), Meru & Bate (2010), and oth- 
ers, find that fragmentation is possible in the outer re- 
gions of protoplanetary disks, but does not occur at the 
5 to 10 AU scale where gas giants are found in the So- 
lar System. Since fragmentation occurs at large radii and 
early times, a large reservoir of mass is typically available 
locally and the likely outcome of fragmentation would 
be very massive planets or brown dwarfs (Stamatellos & 
Whitworth, 2009). 



20 The most important exceptions are the simulations of Boss 
(2008), who finds fragmentation at significantly smaller orbital 
radii than other groups. The origin of this divergent result is not 
clear. 
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3. Comparison with observations 

The architecture of the Solar System's giant planets 
provides qualified support for the core accretion model. 
The time scale for core accretion increases with orbital 
radius, which is qualitatively consistent with the general 
trend of planetary properties in the outer Solar System 

- Jupiter is closest to Solar composition (reflecting a 
fully formed gas giant), while Saturn and the ice giants 
are relatively gas poor. Perhaps these outermost planets 
formed as the gas disk was in the process of being dis- 
persed. Explaining the origin of Uranus and Neptune as 
a consequence of disk fragmentation is not easy. More- 
over the core accretion time scale for the formation of 
Jupiter is reasonable for plausible assumptions. 

Solar System observations also raise doubts about core 
accretion. The time scale to form Neptune, in particular, 
is prohibitively long. This result is now normally inter- 
preted as an indication that Uranus and Neptune may 
not have formed in situ, and as such cannot be used to 
argue against core accretion. It means, however, that the 
ice giants are poor laboratories for testing core accretion. 
Potentially more seriously, measurements of the gravita- 
tional multipole moments of Jupiter (from the Galileo 
orbiter) can be combined with interior structure models 
to yield constraints on the core mass. Until recently there 
was an unambiguous discrepancy between the resulting 
constraints and the predictions of core accretion. Guillot 
(2005), for example, obtained an upper limit on the core 
mass of Jupiter of 10M ffi . For some equations of state 
the constraint on the core mass was below 5M e 21 . This 
was evidently smaller than predictions based on the sim- 
plest models of core accretion (Pollack et al., 1996), and 
in complete agreement with the zero core prediction of 
disk instability. 

This may sound like a clear strike against core accre- 
tion, but in fact matters are not so simple. First, as 
we have already noted fiducial core accretion models are 
based on particular choices of uncertain parameters and 
as such should not be regarded as definitive. Currently, 
it seems reasonable to believe that smaller core masses 

- perhaps as low as 5M© — could be consistent with 
plausible variants of the basic core accretion model (Al- 
ibcrt ct al., 2005b; Hubickyj, Bodenheimer & Lissauer, 
2005; Movshovitz et al., 2010). Second, the "observa- 
tional" constraint on the Jovian core mass is highly in- 
direct. It requires the adoption both of an equation of 
state and of an assumed interior structure model, whose 
parameters are then optimized by comparing the model 
against the actual measured data. There has been sub- 
stantial recent progress in determining the high pressure 
equation of state via ab initio methods, but alas this has 
not eliminated the uncertainty in the quantity we are in- 



terested in - Jupiter's core mass. Recently Militzer et 
al. (2008) have presented Jovian models that include a 
substantial (14 — 18 M©) core, while Nettclmann et al. 
(2008) have computed similarly state-of-the-art models 
that support the Guillot (2005) conclusion that any core 
must be small. The differences appear to result primar- 
ily from different assumptions made by the two groups 
as to the number of distinct layers within the interior of 
Jupiter. NASA's JUNO mission will eventually return 
additional data that will yield new constraints on the 
interior structure of Jupiter, but in the short term there 
appears little prospect of a definitive measurement of the 
planet's core mass. 

Observations of extrasolar planets also yield con- 
straints on the dueling models. Core accretion predicts 
that a greater surface density of planetesimals would lead 
to faster core growth and an increased chance of reaching 
runaway prior to disk dispersal. This is consistent with 
the observed correlation of planet frequency with host 
metallicity (Fischer & Valenti, 2005; Ida & Lin, 2004), 
which appears to be genuinely due to the formation pro- 
cess rather than a metallicity dependence of the migra- 
tion rate (Livio & Pringle, 2003). Taken together with 
the large mass of heavy elements needed to explain the 
small radius of the transiting planet HD 149026 (Sato 
et al., 2005), the observations suggest that most of the 
known exoplanets have properties consistent with the 
outcome of core accretion. 

This does not, of course, mean that disk instability 
does not occur. As we have emphasized, fragmentation 
is expected to occur only at large disk radii, whereas 
almost all known exoplanets have been discovered via 
search techniques that are most sensitive to planets with 
small to intermediate separations. If disk instability does 
occur (and yields planets rather than more massive sub- 
stellar or stellar companions) then we would expect a sec- 
ond population of massive planets in wide orbits (Bolcy 
et al., 2009), with a different host metallicity distribution 
(Rice et al., 2003c). The HR 8799 (Marois et al., 2008) 
system, in which 3 very massive planets orbit at pro- 
jected separations between 24 and 68 AU, is currently the 
only known example of an extrasolar planetary systems 
whose properties hint at a disk instability origin (Dodson- 
Robinson et al., 2009; Krattcr, Murray-Clay & Youdin, 
2010). My own opinion is that the evidence for disk in- 
stability furnished by the HR 8799 system - whose for- 
mation is something of a puzzle for either mechanism - is 
inconclusive. Forthcoming improvements to direct imag- 
ing surveys should yield a much better idea of whether, 
in fact, there is a second population of planets formed 
from fragmentation in the outer disk. 



IV. EVOLUTION OF PLANETARY SYSTEMS 



The same exercise yielded a core mass for Saturn of 1O-2OM0, 
in good accord with the expectations of core accretion 



The story is not over once planets have managed to 
form. Theoretical models, which are now strongly sup- 
ported by observations of the Solar System and of extra- 
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solar planetary systems, suggest at least four mechanisms 
that can lead to substantial post-formation orbital evo- 
lution: 

• Interaction between planets and the gaseous 
protoplanetary disk. This leads to orbital mi- 
gration (Goldreich & Tremaine, 1980) as a con- 
sequence of angular momentum exchange between 
the planet and the gas disk, and can be impor- 
tant for both terrestrial-mass planets and gas giants 
while the gas disk is still present. Gas disk migra- 
tion provides the standard theoretical explanation 
for the existence of hot Jupiters (Lin, Bodenheimer 
& Richardson, 1996). 

• Interaction between planets and a remnant 
planetesimal disk. Planets, especially giant plan- 
ets, can also exchange angular momentum by in- 
teracting with and ejecting planctesimals left over 
from the planet formation process. This mecha- 
nism appears likely to have caused orbital migra- 
tion of at least the ice giants, and possibly also Sat- 
urn, during the early history of the Solar System 
(Levison et al., 2007). 

• Interaction within an initially unstable sys- 
tem of two or more massive planets. There 
is no guarantee that the architecture of a newly 
formed planetary system will be stable over the 
long run. Instabilities can lead to planet-planet 
scattering, which usually results in the ejection of 
the lower mass planets, leaving the survivors on 
eccentric orbits. This could be the origin of the 
typically eccentric orbits seen in extrasolar plane- 
tary systems (Lin & Ida, 1997; Rasio & Ford, 1996; 
Weidenschilling & Marzari, 1996). 

• Tidal interactions between planets and their 
host stars, which are of particular importance for 
hot Jupiters given that their orbital radii are, in 
some cases, just a handful of stellar radii. 

In this Section I discuss each of these mechanisms in turn. 
The focus here is exclusively on dynamical evolution of 
newly formed planetary systems. Of course the surface 
properties of planets also evolve with time, even in the 
absence of orbital evolution, due to changes in the stel- 
lar luminosity and geological processes. Considerations 
of this kind, which are crucial to determining the habit- 
ability of planets, are discussed for example in Kasting, 
Whitmire & Reynolds (1993). 



A. Gas disk migration 

The calculation of the torque experienced by a planet 
embedded within a viscous disk is highly technical, and 
although the basic principles are secure the details are 
anything but. Here, we first give a semi-quantitative 
treatment based on the impulse approximation (Lin & 



Papaloizou, 1979). We then sketch some of the key 
ideas that underly more detailed computations, which 
are based on summing the torque over a set of discrete 
resonances between the planet and the gaseous disk (Gol- 
dreich & Tremaine, 1979). The interested reader is re- 
ferred to Lubow & Ida (2010) for an excellent review of 
the physics of gas disk migration. 



1. Torque in the impulse approximation 

Working in a frame of reference moving with the 
planet, we consider the gravitational interaction between 
the planet and gas flowing past with relative velocity Av 
and impact parameter b. The change to the perpendic- 
ular velocity that occurs during the encounter can be 
computed by summing the force along the unperturbed 
trajectory. It is, 
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where M p is the planet mass. This velocity is directed 
radially, and hence docs not correspond to any angu- 
lar momentum change. The interaction in the two-body 
problem is however conservative, so the increase in the 
perpendicular velocity implies a reduction (in this frame) 
of the parallel component. Equating the kinetic energy of 
the gas particle well before and well after the interaction 
has taken place we have that, 



Av 2 = \6vj_\ 2 + (Av - (5w||) 2 , 
which implies (for small deflection angles), 
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If the planet has a semi-major axis a the implied angular 
momentum change per unit mass of the gas is, 



2G 2 M 2 a 
Aj = b 2 Av 3 ' 
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It is worth pausing at this juncture to fix the sign of the 
angular momentum change experienced by the gas and 
by the planet firmly in our minds. Gas exterior to the 
planet's orbit orbits the star more slowly than the planet, 
and is therefore "overtaken" by the planet. The decrease 
in the parallel component of the relative velocity of the 
gas therefore corresponds to an increase in the angular 
momentum of the gas exterior to the planet. Since the 
gravitational torque must be equal and opposite for the 
planet the sign is such that: 

• Interaction between the planet and gas exterior to 
the orbit increases the angular momentum of the 
gas, and decreases the angular momentum of the 
planet. The planet will tend to migrate inward, 
and the gas will be repelled from the planet. 
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• Interaction with gas interior to the orbit decreases 
the angular momentum of the gas and increases 
that of the planet. The interior gas is also repelled, 
but the planet tends to migrate outward. 

In the common circumstance where there is gas both 
interior and exterior to the orbit of the planet the net 
torque (and sense of migration) will evidently depend 
upon which of the above effects dominates. 

The total torque on the planet due to its gravitational 
interaction with the disk can be estimated by integrating 
the single particle torque over all the gas in the disk. For 
an annulus close to but exterior to the planet, the mass 
in the disk between b and b + db is, 



dm « 2izaYidb 



(215) 



where £ (assumed to be a constant) is some character- 
istic value of the gas surface density. If the gas in the 
annulus has angular velocity f2 and the planet has an- 
gular velocity fl p all of the gas within the annulus will 
encounter the planet in a time interval, 
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Approximating |f2 — f2 p | as 
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which is valid provided that d « s, we obtain the to- 
tal torque on the planet due to its interaction with gas 
outside the orbit by multiplying Aj by dm, dividing by 
Ai, and integrating over impact parameters. Formally 
we would have that, 
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but this integral is clearly divergent - there is what must 
be an unphysically infinite contribution from gas passing 
very close to the planet. Sidestepping this problem for 
now, we replace the lower limit with a minimum impact 
parameter 6 m ; n and integrate. The result is, 



conclude that if all other factors are equal - in particular 
if neither £ in the vicinity of the planet nor 6 m j n vary with 
M p - the time scale for the planet to change its orbital 
angular momentum significantly will scale inversely with 
the planet mass. The migration velocity in this limit will 
then be proportional to M p - more massive planets will 
migrate faster. 

Finally we can estimate the magnitude of the torque for 
parameters appropriate to different stages of giant planet 
formation. First, let us consider a rather low mass core 
(M p = M ffi ) growing at 5 AU in a gas disk of surface 
density £ — 10 2 g cm~ 2 around a Solar mass star. Our 
treatment of the interaction has assumed that the disk 
can be treated as a two-dimensional sheet, so it is ar- 
guably natural to take as a minimum impact parameter 
&min = hfa 0.05a. Using these numbers we find that the 
exterior torque would drive inward migration on a time 
scale, 
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Of course this will be partly offset by the interior torque 
- which has the opposite sign - but absent some physical 
reason why these two torques should cancel to high pre- 
cision we would predict changes to the semi-major axis 
on a time scale of the order of a Myr. This is already 
rapid enough to be a potentially important effect dur- 
ing giant planet formation via core accretion. Second, 
we can evaluate the torque for a fully-formed gas giant. 
A Jupiter mass planet has a Hill sphere m > h, so it 
seems reasonable to argue that the value of 6 m ; n that we 
adopted for an Earth mass core is too small in this case. 
Picking a modestly larger value 6 m i n = 0.2a results in a 
time scale, 



r ~ 2 x 10 5 yr, 



(221) 



that is an order of magnitude shorter than typical pro- 
toplanetary disk lifetimes. If these estimates can be 
trusted to even an order of magnitude the conclusion is 
inescapable - gas disk migration ought to be an impor- 
tant process for the early evolution of the orbits of giant 
planets. 
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(219) 2- Torque at resonances 



It is possible to tidy up this calculation somewhat, for 
example by taking proper account of the rotation of the 
planet frame around the star, and if this is done the result 
is that the expression derived above must be multiplied 
by a correction factor (Papaloizou & Terquem, 2006). 

Aside from the sign of the torque the most important 
result that we can glean from this calculation is that the 
torque on the planet due to its interaction with the gas 
scales as the square of the planet mass. This scaling can 
be compared to the orbital angular momentum of the 
planet, which is of course linear in the planet mass. We 



A more involved - but ultimately more powerful - ap- 
proach to calculate the torque is to decompose it into 
a sum over partial torques exerted at resonant loca- 
tions with the disk (Goldreich & Tremaine, 1979; Tanaka, 
Takeuchi & Ward, 2002). For simplicity, we start by con- 
sidering a planet orbiting a star on a circular orbit with 
angular frequency fl p . A standard exercise in dynamics 
(e.g. Binney & Tremaine 1987) yields the conditions for 
resonances. A corotation resonance exists for radii in the 
disk where the angular frequency fi, 



(222) 
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Lindblad resonances exist when, 

m(VL - Q p ) = ±k (223) 

where m is an integer and kq, the epicyclic frequency, is 
defined as, 

K0 -(^+3O 2 ) (224) 

with $o the stellar gravitational potential. For a Kep- 
lerian potential kq — CI. If we approximate the angular 
velocity of gas in the disk by the Keplerian angular ve- 
locity, the Lindlad resonances are located at, 
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where r p is the planet orbital radius. The locations of 
some of the low order (small to) resonances are shown in 
Figure 28. For an orbiting test particle, the resonances 
are locations where the planet can be a strong perturba- 
tion to the motion. For a gas disk, angular momentum 
exchange between the planet and the gas disk occurs at 
resonant locations. As we noted for the impulse approx- 
imation, the sense of the interaction is that the planet 
gains angular momentum from interacting with the 
gas disk at the interior Lindblad resonances (rj, < r p ). 
This tends to drive the planet outward. The gas loses 
angular momentum, and moves inward. Conversely, the 
planet loses angular momentum from interacting with 
the gas disk at exterior Lindblad resonances (tl, > r p ). 
This tends to drive the planet toward the star. The gas 
gains angular momentum, and moves outward. Notice 
that the gravitational interaction of a planet with a gas 
disk tends — somewhat counter-intuitively — to repel 
gas from the vicinity of the planet's orbit. 

The flux of angular momentum exchanged at each 
Lindblad resonance can be written as, 

T LR (m) ex EM p 2 / c (£) (226) 

where E is the gas density and M p the planet mass. That 
the torque should scale with the square of the planet mass 
is intuitively reasonable — the perturbation to the disk 
surface density scales as the planet mass in the linear 
regime so the torque scales as M 2 . The factor f c (Q is 
the torque cutoff function (Artymowicz, 1993), which en- 
codes the fact that resonances very close to the planet 
contribute little to the net torque. The torque cutoff 
function peaks at, 

^Hm),* 1 (227) 

i.e. at a radial location r ~ r p ± h, where h is the disk 
scale height (this result immediately implies that a three- 
dimensional treatment is necessary for the dominant res- 
onances if the planet is completely embedded within a gas 
disk, as is the case for low mass planets) . The strength of 




FIG. 28 Nominal locations of the corotation (red) and Lind- 
blad resonances (blue) for a planet on a circular orbit. Only 
the low order Lindblad resonances are depicted — there are 
many more closer to the planet. 

the torque exerted at each resonance is essentially inde- 
pendent of properties of the disk such as the disk viscos- 
ity, though of course the viscosity still matters inasmuch 
as it controls the value of the unperturbed disk surface 
density E. 

Figure 29 illustrates the differential torque exerted on 
the disk by the planet, after smoothing over the Lind- 
blad resonances (Ward, 1997). The flux of angular mo- 
mentum is initially deposited in the disk as waves, which 
propagate radially before dissipating. The details of this 
dissipation matter little for the net rate of angular mo- 
mentum exchange. 

Angular momentum transfer at corotation requires ad- 
ditional consideration, and as we will see later thinking 
of these torques in terms of resonances is not as prof- 
itable as for the Lindblad torques. Formally though, in 
a two-dimensional disk the rate of angular momentum 
deposition at corotation is proportional to (Goldreich & 
Tremaine, 1979; Tanaka, Takeuchi & Ward, 2002), 

TaHOcAg) (228) 

where B is the Oort parameter, 

„ r dtt 
= «+- — . (229) 

This implies that in a two-dimensional disk, the reso- 
nant corotation torque vanishes identically in the mod- 
erately interesting case of a disk with a surface density 
profile E cx r~ 3 / 2 . This result does not apply to three- 
dimensional disks (Tanaka, Takeuchi & Ward, 2002). 
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FIG. 29 Schematic illustration of the smoothed torque den- 
sity due to angular momentum exchange between a planet 
and a gas disk at the location of Lindblad resonances, after 
Ward (1997). The peak torque occurs at r ~ r v ±h. The disk 
gains angular momentum from the planet as a result of the in- 
teraction for r > r p , and loses angular momentum for r < r p . 
The interaction is almost invariably asymmetric, such that 
when integrated over the entire disk the planet loses angular 
momentum and migrates inward. 



3. Type I migration 

For low mass planets (generically M p ~ M®, though 
the exact mass depends upon the disk properties) the 
angular momentum flux injected into the disk as a con- 
sequence of the planet-disk interaction is negligible when 
compared to the viscous transport of angular momentum. 
As a result, the gas surface density profile £(r) remains 
approximately unperturbed, gas is present at the loca- 
tion of each of the resonances, and the net torque on the 
planet is obtained by summing up the torque exerted at 
each resonance. Schematically, 

Jplanct = ^ Tlr + 2_j ^ LR TOR (230) 
ILR OLE. 

where the planet gains angular momentum from the inner 
Lindblad resonances (ILR) and loses angular momentum 
to the outer Lindblad resonances (OLR). There is also a 
potentially important co-orbital torque Tcr. Changes to 
the planet's orbit as a result of this net torque are called 
Type I migration (Ward, 1997). 

As noted above (equation 226) the torque exerted at 
each resonance scales as the planet mass squared. If the 
azimuthally averaged surface density profile of the gas 
disk remains unperturbed, it follows that the total torque 
will also scale as M p and the migration time scale, 

r^-^cxM- 1 . (231) 

planet 



Type I migration is therefore most rapid for the largest 
body for which the assumption that the gas disk remains 
unaffected by the planet remains valid. 

Actually evaluating the sum sketched out in equation 
(230) is not easy, and there is no simple physical argu- 
ment that I am aware of that gives even the sign of the net 
torque on the planet. However invariably it is found that 
the Lindblad resonances exterior to the planet are more 
powerful than those interior, so that the net torque due 
to Lindblad resonances leads to inward migration. Note 
that one might think (for example by looking at the sur- 
face density dependence of the torque in equation 226) 
that the sense of migration ought to depend upon the 
surface density gradient — i.e. that a steep surface den- 
sity profile should strengthen the inner resonances rela- 
tive to the outer ones and hence drive outward migration. 
This is not true. Pressure gradients (which depend upon 
the radial dependence of the surface density and tem- 
perature) alter the angular velocity in the disk from its 
Kepler ian value (equation 119), and thereby shift the ra- 
dial location of resonances from their nominal positions. 
A steep surface density profile implies a large pressure 
gradient, so that all the resonances move slightly inward. 
This weakens the inner Lindblad resonance relative to 
the outer ones, in such a way that the final dependence 
on the surface density profile is surprisingly weak (Ward, 
1997). 

Tanaka, Takeuchi & Ward (2002) compute the net 
torque on a planet in a three-dimensional but isothermal 
gas disk. For a disk in which, 

E(r) oc r~ 7 (232) 

they obtain a net torque due to Lindblad resonances only 
of, 

T = -(2.34 - 0.10 7 ) ^D^) a Ep#£ (233) 

This torque would result in migration on a time scale, 
J 

where S p , c s and Q p are respectively the gas surface den- 
sity, gas sound speed, and angular velocity at the loca- 
tion of a planet orbiting at distance r p from a star of 
mass AT*. As expected based on the simple considera- 
tions discussed previously, the migration rate (oc tT ) 
scales linearly with both the planet mass and the local 
disk mass. The time scale becomes shorter for cooler, 
thinner disks — provided that the interaction remains in 
the Type I regime — since for such disks more resonances 
close to the planet contribute to the net torque. 

The most important thing to notice from this formula 
is that the predicted migration time scale is very short. 
If we consider a 5 core growing at 5 AU in a disk 
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with typical parameters (S = 10 2 g cm 2 , h/r = 0.05, 
a = 1) we find, 

tixr ~ 0.5 Myr. (235) 

One concludes that there is a strong argument that 
Type I migration ought to play an important role in the 
formation of giant planet cores. 

The above calculation of the Lindblad torque can be 
verified against hydrodynamic simulations, and is consid- 
ered to be reliable. It is, however, only part of the total 
torque on a planet. What about the co-orbital torque? 
Using a similar linear calculation, it is possible to esti- 
mate the corotation torque as well. This was done by 
Tanaka, Takeuchi & Ward (2002), who derived a total 
torque, 

T = -(1.36 + 0.54 7 ) (^^) 2 (236) 

This formula, and its corresponding migration rate, are 
often described as representing the "standard Type I 
migration rate". One should be aware, however, that 
the formula is derived under conditions (isothermality, 
a smooth density distribution with radius) that do not 
always hold in real disks. It is not, therefore, the com- 
plete answer even in linear theory, and extra caution is 
required before using it under conditions far from its do- 
mian of validity (such as when there is a sharp density 
jump in the disk). 

For very low mass planets, or planets embedded in 
highly viscous disks, the standard Type I migration rate 
calculated from linear theory may be reliable. For higher 
mass planets and / or weaker viscosities, however, it may 
be more profitable (conceptually, and possibly mathe- 
matically) to consider the torque in terms of the dynam- 
ics of closed horseshoe orbits in the co-orbital region. 
These orbits, which are analogous to librating particle 
orbits in the restricted three-body problem, are illus- 
trated in Figure 30. As gas in the disk executes the 
U-shaped turns at the ends of the horseshoe, changes 
in the density of the gas exert a torque on the planet. 
This way of thinking about the torque was considered 
by (Ward, 1991), but largely ignored until simulations 
by Paardekooper & Mellema (2006) uncovered a depen- 
dence of the Type I migration rate on the thermal prop- 
erties of the disk. Subsequently, many authors have stud- 
ied the co-orbital Type I torque in both isothermal (Ca- 
soli & Masset, 2009; Paardekooper & Papaloizou, 2009) 
and non-isothermal (radiative or adiabatic) disks (Klcy, 
Bitsch & Klahr, 2009; Kley & Crida, 2008; Masset & 
Casoli, 2009; Paardekooper et al., 2009). Currently it 
appears as if, 

• The persistence of strong co-orbital torques de- 
pends upon whether or not they are saturated. Sat- 
uration is possible because the region of the disk 
that admits horseshoe orbits is closed and relatively 
small. It cannot absorb or give an arbitrary amount 




FIG. 30 The nonlinear calculation of the torque on an em- 
bedded planet, due to co-orbital gas, is derived from consider- 
ation of the horseshoe drag. The key point is that gas at radii 
close to that of the planet executes closed horseshoe-shaped 
orbits in the corotating frame. Changes in density as parcels 
of gas execute the U-shaped turns result in a non-zero torque 
on the planet. This torque depends on how "non-adiabatic" 
the gas is: does the gas cool radiatively on the time scale 
of the turn? One should also note that the region that sup- 
ports horseshoe orbits is closed. In the absence of viscosity, 
this means that the co-orbital gas can only exchange a finite 
amount of angular momentum with the planet, after which 
the torque saturates. 

of angular momentum to a planet, unless it is "con- 
nected" to the rest of the disk via viscous stresses. 
Large and persistent co-orbital torques are possible 
provided that the disk is viscous enough that the 
torque remains unsaturated. 

• The torque depends upon the cooling time scale 
of the gas in the co-orbital region, measured rel- 
ative to the time required for the gas to execute 
the horseshoe turns. Outward migration under 
the combined influence of co-orbital and Lindblad 
torques (which remain negative) may be possible in 
the inner regions of the disk, where the high optical 
depth results in a long cooling time. 

These results are all very new, and it is reasonable to 
expect that further revisions to our understanding may 
still occur. In particular, little work has yet been done 
to address the question of how realistic turbulent flows 
within the disk affect the torque and its saturation. 

To summarize, Type I migration torques remain poorly 
understood. The co-orbital torque is probably impor- 
tant, but the mathematical relation between the linear 
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FIG. 31 Illustration of the viscous condition for gap opening. 
A gap can open when the time scale for opening a gap of 
width Ar due to tidal torques becomes shorter than the time 
scale on which viscous diffusion can refill the gap. 



calculation, and that done by considering the properties 
of horseshoe orbits, is not clear. It is likely that Type I 
migration is rapid, but the rate and even direction of 
migration may depend upon details of the disk model. 



4. Type II migration 

For sufficiently large planet masses, the angular mo- 
mentum flux from the planet locally dominates the vis- 
cous flux. As a consequence, gas is repelled from high- 



m resonances. The surface density drops near r — 
forming a gap — an annular region in which the surface 
density is smaller than its unperturbed value. 

Two conditions are necessary for gap formation. First, 
the Hill sphere (or Roche radius) of the planet needs to 
be comparable to the thickness of the gas disk, 
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which requires a mass ratio q = M p /M*, 
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FIG. 32 Simulation of the planet-disk interaction in the 
Type II regime in which the planet is sufficiently mas- 
sive to open a gap in the gas disk. Note the presence of 
streams of gas that penetrate the gap region. A movie 
showing the interaction as a function of mass is available at 
http://jilawww.colorado.edu/~pja/planet_migration.html. 



where v = ac s h is the disk viscosity. The time scale to 
open a gap as a result of the tidal torque at an m-th 
order Lindblad resonance is, 



i 2 q 2 fl p 



Ar 



(240) 



Setting t opon = iciosc, and taking m = r p VL p /c s (since, as 
noted above, this value of m is where the torque cutoff 
function peaks), we obtain, 



This condition is satisfied for typical protoplanetary disk 
parameters for q ~ 4 x 10 -4 — i.e. for planet masses 
somewhere between that of Saturn and Jupiter. 

A second condition for gap opening arises due to the 
viscous considerations depicted in Figure 31. To open a 
gap, we require that the tidal torques must be able to 
remove gas from the gap region faster than viscosity can 
fill the gap back in (Goldreich & Tremaine, 1980; Lin & 
Papaloizou, 1980; Papaloizou & Lin, 1984). There are 
various ways to estimate the critical mass above which 
this condition is satisfied. Following Takeuchi, Miyama 
& Lin (1996), we note that the time scale for viscous 
diffusion to close a gap of width Ar is just, 
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For typical disk parameters (h/r = 0.05, a = 10~ 2 ), 
the viscous condition for gap opening is satisfied for q 
modestly larger than 10~ 4 . Combined with the ther- 
mal condition outlined above, we conclude that Jupiter 
mass planets ought to be massive enough to open a gap 
within the disk, whereas Saturn mass planets are close 
to the critical mass required for gap opening. Figure 32 
from Armitage & Rice (2005), shows results from a two- 
dimensional simulation of the planet-disk interaction in 
the Type II regime. Both the gap, and the presence of 
a prominent spiral wave excited within the gas disk, are 
obvious. 
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5. The Type II migration rate 

Once a planet becomes massive enough to open a gap, 
orbital evolution is predicted to occur on the same local 
time scale as the protoplanctary disk. At small orbital 
radii the sense of migration will invariably be inward, 
but the planet will simply follow the motion of the gas 
and can migrate outward in regions where the gas disk is 
expanding (Veras & Armitage, 2004). The radial velocity 
of gas in the disk is, 



M 
' 2irrY, ' 



(242) 



which for a steady disk away from the boundaries can be 
written as, 



v r = — - 
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If the planet enforces a rigid tidal barrier at the outer 
edge of the gap, then evolution of the disk will force the 
orbit to shrink at a rate f p ~ v r , provided that the local 
disk mass exceeds the planet mass, i.e. that irr^T, > 
M p . This implies a nominal Type II migration time scale, 
valid for disk dominated migration only, 
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For h/r = 0.05 and a = 10~ 2 , the migration time scale 
at 5 AU is of the order of 0.5 Myr. 

In practice, the assumption that the local disk mass 
exceeds that of the planet often fails. For example, a 
plausible model of the protoplanetary disk with a mass 
of 0.01 Mq within 30 AU has a surface density profile, 
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The condition that 7rr 2 I] = M p gives an estimate of the 
radius within which disk domination ceases of, 
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Interior to this radius, the planet acts as a slowly moving 
barrier which impedes the inflow of disk gas. Gas piles up 
behind the barrier - increasing the torque - but this pro- 
cess does not continue without limit because the interac- 
tion also deposits angular momentum into the disk, caus- 
ing it to expand (Pringle, 1991). The end result is to slow 
migration compared to the nominal rate quoted above. 
Given a disk model, and in particular a specification of 
how the angular momentum transport efficiency depends 
upon the radius and surface density within the disk, the 
extent of the suppression can be calculated (Ivanov, Pa- 
paloizou & Polnarev, 1999; Syer & Clarke, 1995). To 
illustrate the idea, imagine that we have a disk in which 



the surface density can be written as a power-law in ac- 
cretion rate and radius, 



E oc M a r b , 



(247) 



Syer & Clarke (1995) define a measure of the degree of 
disk dominance, 
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Then for B < 1 (the planet dominated case appropriate 
to small radii) the actual Type II migration rate is (Syer 
& Clarke, 1995), 



T/j = T B- a ^ 1+a \ 



(249) 



Note that with this definition of B, disk dominance ex- 
tends inward a factor of a few further than would be 
predicted based on the simple estimate given above. 

In practice evaluating how the surface density depends 
upon the accretion rate - and thereby determining the a 
which enters into the suppression term - is rather dif- 
ficult, since it requires us to place perhaps unwarranted 
faith in aspects of disk models which are immune from ob- 
servational tests. Proceeding anyway, for the disk mod- 
els of Bell et al. (1997) wc find that a ~ 0.5 at 1 AU 
for M — 10~ 8 Moyr" 1 . At this radius the model with 
a = 10~ 2 has a surface density of about 200 gcm~ 2 . 
For a Jupiter mass planet we then have B = 0.3 and 
Tn = 1-5to. This is a modest suppression of the Type II 
rate, but the effect becomes larger at smaller radii (or for 
more massive planets). It slows the inward drift of mas- 
sive planets, and allows a greater chance for them to be 
stranded at sub-AU scales as the gas disk is dissipated. 
Detailed models suggest that the distribution of massive 
extrasolar planets is consistent with those planets form- 
ing at larger radii, before becoming stranded during the 
migration process due to the dispersal of the gas disk 
(Armitage, 2007; Armitage et al., 2002). 

These estimates of the Type II migration velocity as- 
sume that once a gap has been opened, the planet main- 
tains an impermeable tidal barrier to gas inflow. In fact, 
simulations show that planets are able to accrete gas 
via tidal streams that bridge the gap (Lubow, Siebert 
& Artymowicz, 1999). The effect is particularly pro- 
nounced for planets only just massive enough to open 
a gap in the first place. As the inflowing gas crosses 
the co-orbital region it can exert a substantial additional 
torque on the planet. The Type II migration rate can 
also be qualitatively altered - and even reversed - if two 
planets approach each other in the disk such that their 
gaps start to overlap or such that resonant interactions 
between the planets become important (Masset & Sncll- 
grove, 2001). This effect may have played a role in re- 
ducing the migration of Jupiter and Saturn in the Solar 
System (Morbidelli & Crida, 2007). 
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6. Stochastic migration 

To a first approximation the efficiency of angular mo- 
mentum transport has little impact on the predicted 
Type I migration rate. This assumes, however, that the 
disk is laminar. More realistically, angular momentum 
transport itself derives from turbulence, which is accom- 
panied by a spatially and temporally varying pattern of 
density fluctuations in the protoplanetary disk. These 
fluctuations will exert random torques on planets of any 
mass embedded within the disk, in much the same way 
as transient spiral features in the Galactic disk act to in- 
crease the velocity dispersion of stellar populations (Carl- 
berg & Sellwood, 1985). If we assume that the random 
torques are uncorrelated with the presence of a planet, 
then the random torques' linear scaling with planet mass 
will dominate over the usual Type I torque (scaling as 
Mp) for sufficiently low masses. The turbulence will then 
act to increase the velocity dispersion of collisionless bod- 
ies, or, in the presence of damping, to drive a random 
walk in the semi-major axis of low mass planets. 

To go beyond such generalities, and in particular to es- 
timate the crossover mass between stochastic and Type I 
migration, we need to specify the source of turbulence in 
the protoplanetary disk. MHD disk turbulence driven 
by the magnetorotational instability (Balbus & Haw- 
ley, 1998) provides a well-understood source of outward 
angular momentum transport in sufficiently well-ionized 
disks, and has been used as a model system for study- 
ing stochastic migration by several authors (Laughlin, 
Steinacker & Adams, 2004; Nelson, 2005; Nelson & Pa- 
paloizou, 2004; Yang, Mac Low & Menou, 2009). Density 
fluctuations in MHD disk turbulence have a typical co- 
herence time of approximately half an orbital period, and 
as a consequence are able to exchange angular momen- 
tum with an embedded planet across a range of disk radii 
(not only at narrow resonances). The study by Nelson & 
Papaloizou (2004) was based on both global ideal MHD 
disk simulations, with an aspect ratio of h/r = 0.07, and 
local shearing box calculations. The global runs real- 
ized an effective Shakura-Sunyaev a = 7 x 10~ 3 , which if 
replicated in a real disk would be consistent with obser- 
vational measures of T Tauri disk evolution (Hartmann 
et al., 1998). For all masses considered in the range 
3 < M p < 30 Mg, the instantaneous torque on 
the planet from the MHD turbulent disk exhibited large 
fluctuations in both magnitude and sign. Averaging over 
w 20 orbital periods, the mean torque showed signs of 
converging to the Type I rate, although the rate of con- 
vergence was slow, especially for the lowest mass planets 
in the global runs. These properties are generally in ac- 
cord with other studies of the variability of MHD disk 
turbulence (Hawley, 2001; Winters, Balbus & Hawley, 
2003). Very roughly, the Nelson & Papaloizou (2004) 
simulations suggest that up to M p ~ 10 M s the random 
walk component dominates steady Type I drift over time 
scales that substantially exceed the orbital period. 

How important stochastic (or diffusive) migration is 



for planet formation depends, first and foremost, on the 
strength and nature of the disk turbulence. Many exist- 
ing studies are based on the properties of turbulence sim- 
ulated under ideal MHD conditions, which as we noted 
earlier do not apply to protoplanetary disks. The magni- 
tude of stochastic migration would certainly be reduced 
- though not eliminated - within disks that have a lay- 
ered structure (Oishi, Mac Low & Menou, 2007). Strong 
stochastic migration would pump the mean eccentricity 
(and perhaps inclination) of planetesimals, reducing the 
magnitude of gravitational focusing and potentially lead- 
ing to a greater likelihood of disruptive collisions (Ida, 
Guillot & Morbidelli, 2008). For low mass planets the 
impact of stochastic migration would be to modify their 
survival prospects compared to what would be expected 
from ordinary Type I migration (Johnson, Goodman & 
Menou, 2006). Finally, if stochastic migration still mat- 
ters for rather more massive bodies (M p > M ffi ) then it 
will affect the formation time scale of giant planet cores 
(Rice & Armitage, 2003). 



7. Eccentricity evolution during migration 

Most massive extrasolar planets are observed to be on 
significantly eccentric orbits. Since orbital migration due 
to planet-disk interactions is likely to have occurred in 
these systems, it is of interest to ask whether the same 
process — gravitational interactions between the gas disk 
and an orbiting planet in the Type II regime — also leads 
to excitation of eccentricity. No-one knows for sure. 

The considerations relevant to this problem were set 
out in Goldreich & Tremaine (1980). As with the Type I 
torque, the basic idea is to sum the contributions to e over 
resonances. The number of potentially important reso- 
nances is, however, much larger for an eccentric planet, 
and hence the calculation is (even) harder. Eccentricity 
growth (or decay) depends upon the relative strength of: 

• External Lindblad resonances, which act to excite 
eccentricity. 

• Non-co-orbital corotation resonances, which act to 
damp eccentricity. As noted above, the only coro- 
tation resonance that exists for a planet on a cir- 
cular orbit is co-orbital, so a finite eccentricity is 
necessary for these resonances to be present. 

Unfortunately, the effects leading to damping and exci- 
tation of eccentricity are finely balanced, making robust 
assessment of the sign of the eccentricity evolution dif- 
ficult. The simplest analytic estimates favor damping, 
but only modest saturation of the corotation resonances 
would be needed to tilt the balance in favor of excita- 
tion (Goldreich & Sari, 2003; Masset & Ogilvic, 2004; 
Ogilvie & Lubow, 2003). Numerically, there is general 
agreement that substellar objects of brown dwarf mass 
and above suffer substantial eccentricity growth when 
embedded within a gas disk (Artymowicz et al., 1991; 
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Papaloizou, Nelson & Masset, 2001), while more modest 
excitation has been reported for Jovian mass planets by 
D'Angelo, Lubow & Bate (2006). These results, along 
with related analytic considerations discussed by Moor- 
head & Adams (2008), suggest that although some ec- 
centricity excitation may attend Type II migration it is 
unlikely that the high e tail of the observed exoplanet 
distribution derives solely from this process. 



8. Observational evidence for migration 

The primary interest in Type II migration is as an 
explanation for the existence of hot Jupiters, which are 
observed long after the gas disk has dissipated. If migra- 
tion is common, however, it must be the case that some 
fraction of disk bearing stars harbor embedded migrating 
planets. Unfortunately, no current observational facility 
has the ability to directly image the annular gaps or inner 
holes that are predicted to result from the interaction of 
a planet with the protoplanetary disk, though ALMA has 
a shot at being able to do so (Wolf & D'Angelo, 2005). 
There are, however, a number of T Tauri stars whose 
spectral energy distributions (SEDs) exhibit robust ex- 
cesses in the mid-IR (indicative of gas and dust disks 
at AU scales) without matching excesses in the near-IR 
(Sicilia-Aguilar et al., 2006). Well-known examples of 
such transition disk sources include GM Aur (Calvet et 
al., 2005) and TW Hya (Eisner, Chiang & Hillenbrand, 
2006), but many more such disks have now been identi- 
fied via Spitzer observations (Muzerolle et al., 2010). By 
one common definition, these sources lack optically thick 
inner disks, from which one deduces that small grains are 
absent close to the star, though disks are unquestionably 
present at larger radii. 

What is going on in inner hole sources? Some may gen- 
uinely be stars caught in the act of dispersing their disks 
- perhaps as a result of the photoevaporative mecha- 
nism discussed earlier in these notes. Others, however, 
may be "normal" Classical T Tauri stars around which an 
orbiting planet has created a tidal barrier to the inflow of 
gas and dust, thereby creating an inner hole. Theoretical 
studies suggest that models that invoke the presence of 
planets can fit the observed SEDs (Quillen et al., 2004; 
Rice et al., 2003d), though it is unlikely that this in- 
terpretation is unique or that it can be proved beyond a 
reasonable doubt without spatially resolved observations. 
It is probable that the observed transition sources arise 
from a combination of physical mechanisms, with plan- 
ets being more likely to be implicated in young systems 
that are still accreting gas (Alexander & Armitage, 2007, 
2009). 



their vicinity. The interaction of any remnant planetesi- 
mals with planets, after the dispersal of the gas disk, can 
result in orbital migration of the planets. 

Here, we follow the simple discussion of Malhotra 
(1995) 22 . If we consider a single planetesimal of mass 
Sm interacting with a planet of mass M p at orbital ra- 
dius a there are two possible outcomes, 

• The planetesimal may be scattered outward — pos- 
sibly sufficiently to be ejected — in which case the 
planet moves in by angular momentum conserva- 
tion. Up to numerical factors, 



• The scattering is inward, in which case Sa/a ~ 
+Sm/M p 

Evidently for significant migration to occur we require 
that the total mass in planetesimals be comparable to 
the planet mass, 

^<5m-M p . (251) 

This is a similar result to that obtained in the case of 
gas disk migration, though for planetesimals the restric- 
tion is more severe since while a low mass gas disk can 
drive migration — albeit at a slower pace — ejected plan- 
etesimals are permanently removed from the system and 
cannot influence the planet further. We also note that for 
a single massive planet embedded within a sea of plan- 
etesimals, inward and outward scatterings will at least 
partially balance, leading to little net change in orbital 
radius. 

The foregoing discussion suggests that planetesimal 
migration might be a negligible effect. However, Fernan- 
dez & Ip (1984) showed that the architecture of the outer 
Solar System favors substantial outward migration of the 
ice giants. The key point is that Jupiter is able to eject 
planetesimals from the Solar System more easily that the 
other giant planets. Jupiter itself therefore tends to move 
inward by a relatively small amount due to the ejection 
of debris at initially larger orbital radii. The other outer 
planets scatter bodies inward, to locations from which 
they are removed by Jupiter. This depiction reduces the 
number of outward scatterings, and as a consequence the 
outer planets (minus Jupiter) migrate outward. 

1. Solar System evidence 

Malhotra (1993) and Malhotra (1995) considered the 
effect of the outward migration of Neptune on the ori- 



B. Planetesimal disk migration 

It is unlikely that the formation of gas and ice giant 
planets consumes the entire inventory of planetesimals in 



The treatment here is deliberately over-simplified. The reader in- 
terested in exploring more realistic analytic and numerical mod- 
els is advised to consult Ida et al. (2000) and Kirsh et al. (2009), 
and references therein. 
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FIG. 33 The semimajor axes and eccentricities of known (as 
of 2008) transneptunian bodies. The vertical lines shows the 
location of the 3:2 and 2:1 resonances with Neptune, the 
dashed line shows the minimum eccentricity needed for a body 
to cross Neptune's orbit. 

gin of Pluto and dynamically similar Kuiper Belt Ob- 
jects. The idea is that as Neptune migrated outward, 
Pluto and smaller KBOs were captured into mean mo- 
tion resonances. The eccentricities of captured bodies 
then increase as Neptune continues to move out. For a 
particle locked into a j : j + 1 resonance, the eccentricity 
is (Malhotra, 1995) 

e 2 = e g + — Llnf flNcptunc ) (252) 

.7 + 1 \ a Noptunc,init/ 

where e is the eccentricity on capture into the resonance, 
«Neptunc,init is the semi-major axis of Neptune when the 
particle was captured, and (ZNeptunc is the final semi- 
major axis. For example, if Pluto, then at 33 AU, was 
captured into 3:2 resonance with Neptune when the lat- 
ter was at 25 AU, then migration within the resonance 
out to Neptune's current location at 30.2 AU matches 
Pluto's current eccentricity of e ss 0.25. 

This explanation for the origin of Pluto's peculiar orbit 
is attractive, but even more persuasive evidence for Nep- 
tune's migration comes from the existence of a large pop- 
ulation of KBOs in 3:2 resonance (and smaller numbers 
in other major resonances) with Neptune. This popula- 
tion stands out in even the raw plot of a vs e for KBOs 
shown in Figure 33. In more detail, Murray-Clay & Chi- 
ang (2005) and Hahn & Malhotra (2005) have shown that 
the distribution of KBOs in resonance with Neptune (not 
just the 3:2 resonance) is broadly consistent with, and 
constrains the time scale of, outward migration of Nep- 



tune. Overall, the Solar System evidence seems entirely 
consistent with the hypothesis that substantial migration 
of Neptune captured a substantial disk of planetesimals 
and swept them into resonant configurations akin to that 
of Pluto. 



2. The Nice model 

As noted above, the evidence in favor of Neptune hav- 
ing migrated outward is strong and relatively direct. 
Since the planetesimal scattering that drives this migra- 
tion would also have caused the orbits of Uranus and Sat- 
urn to expand, it is natural to ask whether all the outer 
planets started off in a much more compact orbital config- 
uration. Such a model was proposed by Thommes, Dun- 
can & Levison (1999), who noted that forming Uranus 
and Neptune in the Jupiter-Saturn region of the Solar 
System would immediately alleviate any time scale dif- 
ficulties associated with forming the ice giants in situ. 
Since then, many variations on this theme have suggested 
(Ford & Chiang, 2007; Tsiganis et al., 2005), which differ 
in the details of the planetary initial conditions, the mass 
and structure of the planetesimal disk, and the number 
of planets (possibly more than the current 4) that are 
involved. 

To date, the Nice model (named after the French city) 
of Tsiganis et al. (2005) has received the most attention 
of any specific proposal. The key idea of the Nice model 
is that the early evolution of the outer Solar System in- 
volved two distinct phases, (1) a slow quiescent phase, 
in which planetesimal scattering occurs but the orbits of 
the planets remain almost circular, and (2) a short-lived 
phase of instability in which significant planetary eccen- 
tricities develop. The second phase leads to a brief phase 
of very rapid scattering, which Gomes et al. (2005) as- 
sociate with the Late Heavy Bombardment on the Moon 
about 700 Myr after the formation of the Solar System 
(Hartmann et al., 2000; Strom et al., 2005). Morbidelli 
et al. (2005) show that Jupiter's Trojan asteroids could 
have been captured into their now stable orbits at around 
the same time. 

What dynamics leads to such an evolutionary history? 
In the original version of the Nice model (Tsiganis et 
al., 2005), the idea was that the outer planets started in 
a compact, non-resonant configuration, in which Saturn 
was initially interior to the 2:1 resonance with Jupiter. 
The resonance crossing was then the trigger for the Late 
Heavy Bombardment. Subsequently, Morbidelli et al. 
(2007) have studied a variant in which the planetary 
initial conditions form a resonant chain, so that every 
planet is in a mean-motion resonance with its neighbors. 
In this new version of the Nice model, planetesimal scat- 
tering eventually breaks the resonant chain, and initiates 
the phase of violent instability. "Nice2", as it has been 
dubbed, requires less fine-tuning of initial conditions for 
the planetesimal disk (both variants require a mass in 
planetesimals of about 50 M e ), and appears to provide a 
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better match to the totality of Solar System constraints. 

The Nice model is the closest thing planetary science 
has to a "standard model" that aspires to provide a full 
account of the early history of the outer Solar System. It 
has considerable promise as an economical explanation 
for a diverse array of Solar System observations. As the 
above discussion implies, however, the model-builder has 
considerable freedom to choose initial conditions that are 
broadly consistent with the core accretion paradigm for 
giant planet formation. Balanced against that freedom, 
there are also many potential constraints that can be 
used to rule out models, including, but not limited to, the 
current architecture of the outer and inner Solar System, 
the structure of the Kuiper and asteroid belts, and the 
cratering records on the Moon and other primitive Solar 
System bodies. It is, to say the least, not obvious whether 
these constraints suffice to pin down a unique model. 



y 




-x, x 2 

FIG. 34 Co-ordinate system for the restricted three body 
problem. We work in a co-rotating Cartesian co-ordinate sys- 
tem centered on the center of mass in which the star and 
planet are located at (—xi,0) and (£2,0) respectively. The 
test particle is at position r. 



C. Planet-planet scattering 

While the gas disk is present, gas damping can poten- 
tially protect a multiple planet system against the de- 
velopment of crossing orbits from planet-planet gravita- 
tional interactions (at least if interactions with the gas 
disk actually damp eccentricity, which as noted above is 
somewhat uncertain). Once the gas is gone, gravity can 
go to work on what may be an unstable planetary sys- 
tem and change the orbital radii and eccentricities of the 
planets. This process — gravitational scattering — is 
probably the most widely invoked mechanism to explain 
the large eccentricities of many extrasolar giant planets. 



1. Hill stability 

Let us begin with some analytic considerations. The 
general N-body problem of the motion of N point masses 
interacting under Newtonian gravity is analytically insol- 
uble for N > 2. Here, we start by considering a special 
case of ./V = 3 in which two bodies, of arbitrary mass, 
have a circular orbit, while a third body of negligible mass 
orbits in the known gravitational field of the massive ob- 
jects. This problem — called the circular restricted 3- 
body problem — still defies analytic solution, but it is 
possible to place useful limits on the motion of the third 
body (often described as a "test particle"). The circular 
restricted 3-body problem is a reasonable approximation 
to several situations of great practical interest, including 
the motion of asteroids in the vicinity of Jupiter, and 
the evolution of planetesimals near a growing planet. A 
good description of the problem can be found in Murray 
& Dermott (1999), whose treatment we largely mirror 
here. The more general 3-body problem is discussed (in 
both the planetary and multiple star contexts) in a book 
by Valtonen & Karttunen (2006). 

As shown in Figure 34, we consider a binary system 
in which the massive bodies have mass m\ and mi re- 



spectively. We work in a corotating co-ordinate system 
centered on the center of mass. The orbital plane is (x, y) 
in Cartesian co-ordinates, and the test particle is located 
at position r. 

If the angular velocity of the binary is f2, the equations 
of motion for the test particle are, 

r = -V$ - 2 (ft x r) - ft x (ft x r) (253) 

$ = -^l-^i. (254) 
n r 2 

Expressed in components, we have, 
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The acceleration due to the centrifugal force can be sub- 
sumed into a pseudo-potential. Defining, 
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Digressing briefly, we note that U is (up to an arbitrary 
minus sign) the "Roche potential" . Two stars, or a star 
plus a planet, that rotate synchronously while on circular 
orbits occupy Roche equipotentials. If their size is com- 
parable to the size of the Roche lobe — defined by the 
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critical figure-of-eight shaped equipotential that passes 
through the inner Lagrange point Li — then the bodies 
suffer significant tidal distortion. A useful approximation 
for the radius -Rrl of a sphere with the same volume as 
the Roche lobe was provided by Eggleton (1983). For a 
binary with mass ratio q = majrax and separation a, 



0.49g 



2/3 



Rrl ^ 

~a~ ~ a6(7 2 / 3 + ln(l + qV 3 ) 



(258) 



This equation can be used to assess, for example, how 
close hot Jupiters are to overflowing their Roche lobes. 
For a Jupiter mass planet with q = 10 -3 , 



i? RL ~ 0.048a. 



(259) 



A planet with the same radius as Jupiter (7.14 x 10 9 cm) 
would then overflow its Roche lobe interior to a = 
0.01 AU. A very short period hot Jupiter, such as OGLE- 
TR-56b (Torres et al., 2004) with a period of 1.2 days, 
has a semi-major axis that is about 0.0225 AU. So this 
planet, and more securely other hot Jupiters that or- 
bit modestly further out, is safe against mass transfer, 
though not by a large margin. 

Returning to the general equations (257), we eliminate 
the Coriolis terms by multiplying through by x, y and z 
and adding. We then obtain, 
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(260) 



where v is the velocity and Cj, called the Jacobi constant, 
is the arbitrary constant of integration. Cj is an energy- 
like quantity that is a conserved quantity in the circular 
restricted 3-body problem. 

The existence of this integral of motion is important 
because it places limits on the range of motion possible 
for the test particle. For a particle with a given initial po- 
sition and velocity, we can use equation (260) to compute 
Cj, and hence to specify zero-velocity surfaces, defined 
via, 



2U = Cj, 



(261) 



which the particle can never cross. If the volume enclosed 
by one of the zero- velocity surface is finite, then a particle 
initially within that region is guaranteed to remain there 
for all time. This concept is known as Hill stability. 

The topology of the zero-velocity surfaces in the re- 
stricted three-body problem varies according to the value 
of Cj. An example is shown in Figure 35. In this instance 
the zero-velocity surfaces define three disjoint regions in 
the (x, y) plane, one corresponding to orbits around the 
star, one corresponding to orbits around the planet, and 




FIG. 35 Forbidden zones (dark regions) in an example of 
the restricted 3-body problem. For this particular choice of 
the Jacobi constant Cj, particles can orbit the star at small 
radii; the planet in a tight orbit; or the star-planet binary 
as a whole. The existence of zero-velocity surfaces, however, 
means that particles cannot be exchanged between these re- 
gions. 



one corresponding to orbits around the star-planet bi- 
nary. A particle in any one of these states is stuck there 
— it cannot cross the forbidden zone between the differ- 
ent regions to move into a different state. 



2. Scattering and exoplanet eccentricities 

The test particle analysis discussed above can, some- 
what surprisingly, be extended to the much tougher prob- 
lem of the stability of two planets orbiting a star. Con- 
sider the situation shown in Figure 36, in which planets 
of mass ?7i2 and 7713 orbit a star of mass mi in circular 
orbits with semi-major axes 02 and 03 respectively. The 
stability of the system evidently must depend upon the 
relative, rather than the absolute, spacing between the 
orbits. Accordingly we write, 



a 3 = a 2 (l + A) 



(262) 



with A being a dimensionless measure of the orbital sep- 
aration between the planets. We further define /_t2 — 
m2/mi and 113 = 7713/7774. Then for ^2,^3 ^ L Glad- 
man (1993), drawing on earlier results derived by Mar- 
chal & Bozis (1982) and others, showed that the system 
is guaranteed to be stable provided that the separation 
A exceeds a critical separation A c given by, 



2.40 (M2 + M3) 



1/3 



(263) 
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FIG. 36 Setup for the stability calculation of a two planet 
system in which both of the planets are on circular orbits. 
Unlike in the case of the Hill problem, here we strictly require 
that mi -C mi and 7713 <C mi. 

Note that analytic results leave open the question of 
whether systems with A < A c are actually unstable, all 
we know is that A > A c is sufficient for stability. This 
condition reduces to the test particle result if ^3 — > 0, 
as of course it should 23 . As an example, if we compute 
the critical separation for planets of the mass of Jupiter 
and Saturn, we obtain A c ~ 0.26. The actual separation 
of Jupiter and Saturn in these units is A ~ 0.83, so an 
isolated planetary system in which Jupiter and Saturn 
were on circular orbits would assuredly be stable for all 
time. 

What about more complex systems? It is possible to 
include non-zero eccentricities into this analysis, but not 
more planets. For a multiple planet system one might 
plausibly reason that the system will be unstable if any 
pair substantially violates the critical two-planet separa- 
tion for Hill stability. It is also true that the system will 
generally become more stable as the separations increase 
(Chambers, Wetherill & Boss, 1996). However, no abso- 
lute stability bound is known for any planetary system 
with N > 3. 

If a two-planet system is unstable, the possible out- 
comes of the instability can be divided into four classes: 

1 . The separation evolves (increases) until the system 



Note, however, that the analysis for the restricted three-body 
problem applies for an arbitrary mass ratio of the massive bodies, 
whereas the result for two planets requires that both be much less 
massive than the star. 



achieves a state that is stable over the long term. 

2. One planet is ejected, while the other remains 
bound, generally with e ^ 0. 

3. The planets physically collide. 

4. One planet impacts the star, or is scattered into a 
short-period orbit for which tidal effects are impor- 
tant. 

The last two channels are not possible in a model 3-body 
problem, in which the planets are represented by point 
masses, but can occur (especially planet-planet collisions, 
which become frequent at small radii) in real systems. 

The idea that gravitational scattering and planetary 
ejections might account for the eccentricity of extraso- 
lar planets was proposed as soon as it became clear that 
extrasolar planets were not typically on circular orbits 
(Lin & Ida, 1997; Rasio & Ford, 1996; Weidenschilling & 
Marzari, 1996). Quantitative study of such models re- 
quires large-scale N-body integrations, first to derive the 
statistical distribution of outcomes of any given scenario 
(since the systems are typically chaotic, nothing can be 
said about any single run), and second to map out the 
large parameter space that results when one considers 
different numbers of planets with different initial separa- 
tions, masses and so forth. 

Ford, Havlickova & Rasio (2001) presented a compre- 
hensive study of the dynamics of equal mass two planet 
systems. The planets were set up on circular orbits close 
to the stability boundary, and allowed to evolve under 
purely N-body forces until the system relaxed to a stable 
state. They found that the predicted fraction of collisions 
increases sharply for small orbital radii and / or larger 
planetary radii. For pairs of Jupiter mass and Jupiter 
radius planets initially located at 5 AU, the most com- 
mon outcome is two planets (65%), followed by ejections 
(35%), with collisions (10%) a distant third. If the same 
pair of planets starts at 1 AU, however, collisions occur 
roughly 30% of the time. This conclusion is important 
for studies of extrasolar planet eccentricity, because col- 
lisions yield relatively low eccentricities for the merged 
planet. Indeed, Ford, Havlickova & Rasio (2001) found 
that equal mass planet scattering failed to match the 
observed distribution of eccentricities. However, subse- 
quent calculations that relaxed the equal mass assump- 
tion showed that two planet systems in which the planets 
have a realistic range of masses can yield agreement with 
observations (Ford, Rasio & Yu, 2003). 

There is only a rather small range of orbital separations 
which allows a two planet system to be unstable over the 
long term (greater than around 10° yr, which is roughly 
the dispersal time for the gas disk), while not being vio- 
lently unstable. This observation means that it is easier 
to set up an internally self-consistent scattering model 
with three or more planets, since a wider range of such 
systems eventually lead to interesting dynamics. Models 
starting with three or more planets have also been studied 
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FIG. 37 The differential (upper panel) and cumulative (lower 
panel) eccentricity distribution of known extrasolar plan- 
ets (grey curves) is compared to the predicted distribution 
that results from scattering in three planet systems (shown 
in green). The simulation results were derived by numer- 
ically evolving an ensemble of unstable planetary systems 
made up of three planets whose masses were drawn from the 
observed mass function for extrasolar planets in the range 
0.3 M.j < M p < 5 Mj. The inner planet was initially at 
a — 4.5 AU. Based on simulations by Raymond et al. (2008), 
updated with current data as of Spring 2010. 



in some detail (Adams & Laughlin, 2003; Marzari & Wei- 
denschilling, 2002; Terquemfe Papaloizou, 2002). Recent 
comprehensive studies, such as those by Chatterjee et al. 
(2008) and Juric & Tremaine (2008), find that scattering 
models yield a good quantitative match to the observed 
distribution of extrasolar planet eccentricities. Just how 
good this agreement is is illustrated in Figure 37, which 
shows how the final simulated eccentricities compare to 
the data (Raymond et al., 2008). Although this match 
does not prove that planet-planet scattering is the sole 
(or even the main) source of exoplanet eccentricities, it 
is sufficient to establish this model as the leading candi- 
date for explaining the broad distribution of exoplanet 
eccentricity. 



D. Predictions of migration theories 

In summary, there is persuasive circumstantial evi- 
dence for the action of at least three separate processes 
that lead to the early evolution of planetary systems: 

• Gas disk migration in the Type II regime ap- 
pears to be necessary to explain the existence of 
hot Jupiters. An argument can be made that more 



modest levels of migration ought to be common, 
since at least some gas must necessarily be present 
at the epoch when giant planets form. 

• Planetesimal disk migration provides a persua- 
sive explanation for the origin of Pluto's odd orbit 
together with some of the more detailed proper- 
ties of the Kuiper Belt. One can make a case that 
this process too ought to be common in the outer 
reached of planetary systems. Gas giant formation 
almost certainly becomes more difficult further out 
in the disk, so it is quite plausible that the zone 
where gas and ice giants manage to form is often 
surrounded by a disk of planetesimal debris that 
has been unable to grow to large sizes. 

• Planet-planet scattering works well as an ex- 
planation for the eccentricity distribution of giant 
extrasolar planets. There is no straightforward in- 
dependent argument that the unstable initial con- 
ditions needed for such models to work are generi- 
cally realized in nature, but the empirical evidence 
seems to suggest that they are. 

A number of qualitatively different tests of these the- 
oretical ideas are possible in the near future. Planet- 
planet scattering models, for example, predict that the 
surviving planets ought to have a distribution of incli- 
nation as well as eccentricity, which can be measured 
(relative to the stellar equatorial plane) for those objects 
observed in transit (Rossiter, 1924; Winn et al., 2009). 
With a large sample of measured inclinations it ought 
to be possible to determine whether the predictions of 
pure planet-planet scattering models are confirmed, or 
whether the final inclinations are instead also affected by 
other processes, such as damping from a residual gas disk 
or Kozai pumping of inclination by more distant com- 
panions (Nagasawa, Ida & Bessho, 2008). Planet-planet 
scattering also predicts the existence of a (small) pop- 
ulation of very weakly bound, typically eccentric plan- 
ets, with semi-major axes of 10 2 AU and more (Scharf & 
Menou, 2009; Veras, Crepp & Ford, 2009). Direct imag- 
ing surveys of young stars have the potential to detect 
this distinctive population. 

The combined action of multiple evolutionary mecha- 
nisms may also give rise to new classes of planetary sys- 
tems. At larger orbital radii (than those currently probed 
by observations of exoplanets) it seems likely that we 
ought to see planetary systems whose dynamics has been 
affected by both planet-planet scattering and planetesi- 
mal disk migration. N-body simulations suggest that the 
signature of this combination is a transition from gener- 
ally eccentric to nearly circular planetary orbits as the 
mass of the planetary system is reduced (Raymond, Ar- 
mitage & Gorelick, 2009, 2010). If true, the near-circular 
orbits of the giant planets in the Solar System might in 
fact be typical of the architecture of relatively low-mass 
systems at large orbital radii. For higher mass systems 
the same simulations predict a high abundance of reso- 
nant configurations, including resonant chains that would 
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be planetary analogs of the Laplace resonance in the Jo- 
vian satellite system. 



E. Tidal evolution 

Two bodies in a tight orbit experience tidal forces as 
a consequence of the gradient in the gravitational force 
across their finite radius. The tidal forces raise tidal 
bulges on the surface of the bodies, whose shape is ap- 
proximately defined by the condition of hydrostatic equi- 
librium in the asymmetric gravitational potential. If the 
axis of the tidal bulge is offset with respect to the line 
joining the centers of the two bodies, the result is a tidal 
torque which modifies the semi-major axis and eccentric- 
ity of the system. 

Tides are dynamically important in the Solar System. 
Energy dissipation associated with oceanic tides, raised 
on the Earth by the Moon, is responsible for a slow 
but measureable increase in the Earth-Moon separation 
(Dickey et al., 1994). The basic framework for under- 
standing these phenomena dates back to work by George 
Darwin (Charles' son) more than a century ago (Darwin, 
1879), but a general first-principles theory of tides re- 
mains elusive. As a result, although it is safe to say that 
tidal effects must be important for hot Jupiters (and have 
very probably been responsible for circularizing their or- 
bits), the quantitative evolution of hot Jupiters due to 
tides remains subject to substantial uncertainty. 




FIG. 38 Upper panel: illustration of the hydrostatic tidal re- 
sponse of a body (a star or a planet) of mass M and radius R 
to a point mass companion of mass m orbiting at distance a. 
A tidal bulge of amplitude £ is raised on the surface, aligned 
with the separation vector between the two bodies. Lower 
panel: in the presence of dissipation, the tidal bulge lags (or 
leads, depending upon the spin) the motion of the compan- 
ion by an angle cj>. As a consequence, the tidal bulge - now 
idealized as two point masses of mass AM - exerts a torque 
which modifies the orbit of the secondary. 



1. The tidal bulge and tidal torque 

Even the elementary theory of tides is quite intricate. 
We can gain considerable physical insight, however, from 
a simple "back of the envelope" calculation that ignores 
innumerable order unity numerical factors and effects. 
Consider the tide raised on a fluid body of mass M and 
radius R by a companion, of mass m, orbiting in a cir- 
cular orbit at distance a. The geometry is shown in 
Figure 38. We seek to determine, first, the height of 
the tide £, and, second, the torque that results if the 
tidal bulge is misalinged by some angle cf> with respect to 
the line joining the centers of the two bodies. We will 
assume, throughout, that the tidal deformation corre- 
sponds closely to the hydrostatic response of the body 
in the gravitational field defined by both bodies (the 
"static" tide). This is a reasonable approximation. Ap- 
plying the virial theorem to the body on which the tide 
is raised, we find that the central sound speed ought to 
be comparable to the orbital velocity around the body 
at radius R. For a ~> R, it follows that one orbit of the 
companion corresponds to many sound crossing times of 
the fluid, and hydrostatic equilibrium has time to be es- 
tablished. 

To estimate the height of the tidal bulge, £, we note 
that the gravitational force (per unit mass) exerted by 
the companion on the near side of the body differs from 



that exerted at the center by an amount, 
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This tidal force allows us to raise a tidal bulge up to a 
height where the self-gravity of the fluid body is reduced 
by the same amount, 



d (GM\ Gm 



dR V R 2 J 



R 



i 
R 



/m\ ( R 
\MJ la 



(265) 
(266) 



The height of the bulge falls off rapidly for a ^> R. If 
we assume that the fluid body has a uniform density 
p (and thereby ignore important structural factors) the 
mass associated with the bulge is given by, 



AM = 4ttR 2 £ P , 
which simplifies to yield, 
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The work required to raise the bulge against the self- 
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gravity of the fluid body is, 
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This is the amount of energy associated with the tidal 
deformation of the body. 

Let us now backtrack to revisit the assumption that the 
tidal response of the fluid is hydrostatic. If this were ex- 
actly true, the tidal bulge would line up precisely along 
the line joining the centers of the two bodies, and, by 
symmetry, there would be no tidal torque (as shown in 
the upper panel of Figure 38). In reality, however, the 
tide represents the response of the fluid to forcing at some 
non-zero frequency f2, given (for a companion on a circu- 
lar orbit around a non-rotating primary) simply by the 
orbital frequency. If the central body is non-rotating (or 
rotating slower than the orbital frequency), then the de- 
parture from hydrostatic equilibrium occasioned by the 
finite response time will cause the bulge to lag by some 
angle (p. The bulge will lead if the spin of the central 
body is faster than the orbital frequency. A lagging or 
leading bulge will exert a torque on the companion that 
cause the orbit to decay or expand in radius, respectively 
(e.g. the Moon's orbital period exceeds the Earth's spin 
period, so the torque for the Earth-Moon system results 
in recession of the lunar orbit away from the Earth). 

We now assume that the torque due to the tide can be 
represented as that due to two point masses AM offset 
from the line of centers by a small angle <j> (Figure 38, 
lower diagram). The tidal force is then, roughly, 
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and the work done by the tidal force per orbit of the 
companion is, 
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If we knew - or could calculate - <f>, this equation (with 
the missing numerical factors restored) would yield the 
rate of decay or rate of expansion of the orbit due to the 
tides. 

Calculating (j> from first principles is a hard task. We 
can, however, at least gain a more transparent under- 
standing of the physical problem of what determines <j). 
To do so, we define the tidal Q as the ratio of the en- 
ergy stored in the tidal deformation to the energy that is 
dissipated in one cycle, 
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Using our estimates, we have that, 
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This is an important result. We see that the magnitude of 
the tidal lag, and hence the strength of the tidal torque, 
is directly linked to the amount of dissipation within the 
tidally distorted body. 

Finally we can estimate the rate of decay of the orbital 
separation due to the (lagging) tidal bulge. For m <C M 
we have that, 
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where P = 2-Ky / a 3 /(GM) is the orbital period. Substi- 
tuting for the various quantities, we obtain, 
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Up to factors of the order of unity, this expression agrees 
with the standard tidal formula quoted, for example, by 
Jackson, Greenberg & Barnes (2008), who also give the 
corresponding expressions for the change in orbital eccen- 
tricity. Standard references for the astrophysical theory 
of tides include Goldreich & Soter (1966) and Hut (1981). 



2. Determining the tidal Q 

The above analysis suffices to illustrate two important 
points, 

• Tidal forces decline extremely rapidly with increas- 
ing orbital separation. 

• The rate of tidal evolution depends upon the 
amount of dissipation present within the two bodies 
that are interacting tidally. 

The attentive reader may also have noticed the astro- 
physical sleight of hand by which all manner of in- 
tractable physics has been swept into a single unknown 
parameter, Q. Approaches to estimating Q can be di- 
vided into those that rely on extrapolations from mea- 
sured values in well-observed systems, and those that at- 
tempt to compute Q by identifying dissipative processes 
within stars of planets. 

For giant planets, the primary observational constraint 
on Q comes from direct measurements of the tidal evo- 
lution in the orbits of the Galilean satellites (Lainey et 
al., 2009; Yoder & Peale, 1981) 24 . These empirical es- 
timates suggest that Q ~ 10 5 , and based on this may 



24 Historically, estimates of Jupiter's Q have primarily been derived 
indirectly, by assuming that the excess heat flux from Io derives 
from tidal effects. Greenberg (2010) gives a good review of these 
estimates. 
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workers adopt values for the tidal Q of extrasolar planets 
that are similar (typically in the range Q = 10 5 — 10 6 ). 
Considerable caution as to the validity of this extrapo- 
lation from Jupiter to extrasolar planets is, however, in 
order. As noted already, the Q of a star or planet is not 
some fixed and immutable property of a body akin, say, 
to its mass. Rather, Q describes the response of a body 
to forcing at one or more specific frequencies (in the case 
of the Jovian estimate, the frequency of relevance is the 
difference between the spin frequency of Jupiter and the 
orbital frequency of Io). To extrapolate correctly to ex- 
trasolar planets, we need to make some assumption as to 
how Q varies as a function of frequency. One simple as- 
sumption is to postulate that the tidal lag angle remains 
constant - in which case Q is indeed a constant - but one 
might equally assume that the tidal bulge has a constant 
time lag, in which case Q — Q(f2). If one is interested in 
tidal eccentricity evolution, moreover, the forcing is not 
monochromatic but rather has a spread across a range 
of frequencies. As a result of these complications, em- 
pirical models of extrasolar tidal evolution are subject to 
substantial but unquantifiable uncertainties. 

Given the uncertainties in the empirical approach, a 
theoretical determination of Q would evidently be ex- 
tremely valuable. Achieving this goal requires first iden- 
tifying, and then calculating, the primary source of dissi- 
pation that acts on the tide. Molecular viscosity is insuf- 
ficent, so we are left with a variety of hard-to-calculate 
candidates that include non-linear dissipation of waves 
and interactions between the tide and turbulent processes 
within the body. Recent theoretical work (Goodman & 
Lackner, 2009; Ogilvie & Lin, 2004) in this area, although 
it still falls short of being able to predict Q from first- 
principles, has nonetheless proven influential in identify- 
ing additional properties of planets that may influence 
the Q. The presence of a rigid, solid core, for example, 
can substantially alter the tidal response of a planet as 
compared to an observationally almost indistinguishable 
body lacking a core. Given the rapidly improving obser- 
vations of extrasolar planets that are surely vulnerable 
to tidal evolution, one may hope that this is an area ripe 
for further theoretical and observational progress. 
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